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ABSTRACT OF THE THESIS 


ANALYTICAL STUDY 
of the 

LIQUID PHASE TRANSIENT BEHAVIOR 

of a 

HIGH TEMPERATURE HEAT PIPE 
by 

Gregory Lawrence Roche 
Master of Science in Mechanical Engineering 
University of California, Los Angeles, 1988 
Professor Ivan Catton, Chair 

The transient operation of the liquid phase of a high temperature heat pipe 
is studied. The study was conducted in support of advanced heat pipe 
applications that require reliable transport of high temperature, high flux 
thermal energy over small temperature drops and significant distances under 
a broad spectrum of operating conditions. 

The heat pipe configuration studied consists of a sealed cylindrical enclosure 
containing a capillary wick structure and sodium working fluid. The wick is 


an annular flow channel configuration formed between the enclosure interior 
wall and a concentric cylindrical tube of fine pore screen. 

The study approach is analytical through the solution of the governing 
equations. The energy equation is solved over the pipe wall and liquid region 
using the finite difference Peaceman-Rachford Alternating Direction Implicit 
numerical method. The continuity and momentum equations are solved over 
the liquid region by the integral method. The energy equation and liquid 
dynamics equations arc tightly coupled due to the phase change process at the 
liquid-vapor interface. A kinetic theory model is used to define the phase 
change process in terms of the temperature jump between the liquid-vapor 
surface and the bulk vapor. Extensive auxiliary relations, including sodium 
properties as functions of temperature, are used to close the analytical system. 

The solution procedure is implemented in a Fortran computer algorithm 
with some optimization features to take advantage of the IBM System/370 
Model 3090 vectorization facility. Code output includes: temperatures of the 
pipe wall and liquid; velocities and pressures of the liquid; radii of curvature 
of the liquid-vapor interface meniscii; input, output, evaporation and 
condensation heat fluxes; sonic, entrainment, and capillary heat flux limits; 
and identification of operating limit violations. 

The code was intended for coupling to a vapor phase algorithm so that the 
entire heat pipe problem could be solved. As a test of code capabilities, the 
vapor phase was approximated in a simple manner. The vapor was assumed 
to be quasisteady relative to liquid dynamics and uniform throughout specified 
control volumes. The simple treatment allowed demonstration of calculations, 


xvn 



but did not necessarily produce physically significant results. The kinetic 
theory phase change model was found to produce numerical stability and 
accuracy problems. The solution approach would also benefit from special 
treatment of the widely disparate length scales in the different coordinate 
directions. Cross-stream temperature gradients were found to be small such 
that a simplified solution approach taking advantage of this feature could be 
developed. 
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Chapter I 

INTRODUCTION 


a heat pipe may be regarded as a synergistic engineering structure 
which is equivalent to a material having a thermal conductivity greatly 
exceeding that of any known metal.''' 

G. M. Grover et al. 1964 [/] 


1.1 THE HEAT PIPE 

Transporting thermal energy at a high rate over a small temperature 
gradient is an important requirement for many heat transfer applications. One 
common method for meeting this requirement involves evaporation of a liquid 
at the heat source, transporting the resulting vapor to a heat sink for 
condensation, and replenishing the liquid through return of the condensate or 
provision of a continuous fresh liquid supply. The disadvantage of this 
approach is that generally work in the form of gravity or mechanical pumping 
must be performed on the system to maintain the supply of liquid. The heat 
pipe provides a passive alternative heat transfer system that accomplishes the 
same heat transfer functions without requiring external work. 
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The heat pipe generally consists of a long sealed container that encloses a 
working fluid and a capillary wicking structure as shown in Figure 1. Heat 
added at the evaporator section causes a temperature increase and evaporation 
of the working fluid liquid phase. Mass addition to the vapor space and the 
temperature rise produce a vapor pressure gradient resulting in vapor flow 
from the heated evaporator to the cool condenser where condensation occurs. 
The condensate is returned to the heated end by the action of capillary forces 
in the wicking material. Since the primary mode of heat transfer is through 
the latent heat of vaporization, steady state operation is nearly isothermal even 
while maintaining a high heat transfer rate. 
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1.2 THE PROBLEM 


Heat pipes are self-starting, self-regulating, high efficiency devices that 
operate passively. Unlike normal conducting materials, however, a heat pipe 
cannot be characterized by a single property such as an effective thermal 
conductivity. Rather, operation of the heat pipe is governed by the coupled 
laws of fluid dynamics and heat transfer as applied to each unique 
configuration. 

There have been many advances in heat pipe design, theory and practice 
since the first experimental work performed by Grover in the early 1960's and 
the first quantitative analysis performed by Cotter in 1965[1][2] . There are 
now over 2000 heat pipe references in the literature[3] . There have been six 
international heat pipe conferences, and there are at least five heat pipe 
textbooks and design guides[4][5][6][7][8] . The scope of published 

information spans a broad range of theory, analysis, experiment, and 
application. However, the majority of these works are limited to steady state 
applications under highly idealized conditions. Comparatively little work has 
been done that addresses future heat transport requirements for applications 
that would benefit from the use of heat pipes. For example, space missions 
requiring high power supplies and reliable operation over a wide range of 
conditions may use nuclear generating systems. Such a system has a wide 
range of heat transfer requirements as shown in Figure 2 that may be met by 
heat pipes[9] . The heat transfer requirements range from high temperature, 
high flux primary heat transport to lower temperature secondary transport to 
heat transport for support equipment. 
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Figure 2: Heat Pipe Applications in a Space Nuclear Power System 

Adapted from Merrigan [9]. 


Primary heat transport channels reactor thermal energy to the power 
conversion system for production of electricity. Waste heat rejection transfers 
unuseable thermal energy to a radiator system for disposal to space. Support 
equipment heat transfer requirements include electromagnetic pump cooling, 
radiation shield cooling, control drum and motor cooling, and startup and 
shutdown thermal conditioning of liquid metal pumped heat transfer loops. 
As detailed designs are developed heat pipes may also be used in the control 
of power conditioning equipment, in control electronics cooling, and in low 
temperature waste heat rejection. 

Heat transfer requirements unrelated to space nuclear power that may be 
met by heat pipes have also been identified. The leading edges of hypersonic 
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re-entry vehicles could be cooled and maintained isothermal by liquid metal 
heat pipes[10] . While the working fluid would normally be below the melting 
point, re-entry heat loads would melt the working fluid and initiate heat pipe 
operation. High temperature heat pipes have also been proposed for cooling 
and making isothermal the nozzles of rocket engines[l 1] . 

Heat pipes may be the preferred design approach for many of these 
applications. Heat pipes offer a highly redundant system of independent flow 
paths. Passive operation can provide automatic handling of reactor decay heat 
as well as for restart after extended periods of inoperation without the need for 
stored power for pumps and controls. Heat pipes operate with minimal 
temperature drop such that high efficiency thermal radiation can be achieved. 
Heat pipes are able to transform thermal energy between high flux and low 
flux conditions at the pipe external boundaries. 

While the range of benefits is impressive, many issues require further 
research. One of the most challenging issues is the transient operation of heat 
pipes. Transient heat pipe behavior has become the focus of recent research 
due to requirements for reliable operation over wide operating conditions. 
Consider the space-based nuclear power generating system. The system would 
be placed in orbit at temperatures below the melting point of the heat pipe 
working fluid. Placing the system in operation would require heat pipe startup 
from the frozen working fluid state. After achieving steady state operation, 
nominal changes in reactor power production, power conversion system 
operation, or heat rejection efficiency would result in transients requiring heat 
pipe response to load level changes. At some point due to maintenance 
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requirements or duty cycle, system shutdown would be required. The 
shutdown may be to temperatures below the working fluid melting point. The 
heat pipe must attain a state that permits eventual restart from the frozen or 
liquid state. The heat pipe may be required to operate for extended periods in 
a standby or low power level mode with occasional short bursts of full power 
operation. 

Issues that must be resolved include determining the heat pipe response 
time, the ultimate capability of the heat pipe to handle the transient, and 
whether future operation is possible if the heat pipe is taken outside of design 
capabilities. The useable lengths of heat pipes must be determined, as well as 
optimal design configuration. Durability is an issue due to the cyclic nature 
of the application. Finally, the heat pipe problem must be reduced to a 
minimal set of parameters suitable for a system study with high fidelity results. 


1.3 THE OBJECTIVES 

The purpose of this research is to develop an analytical model of the 
transient dynamics of the liquid phase in a high temperature heat pipe. This 
liquid phase model is intended to be coupled with a vapor phase model for the 
assessment of the potential for heat pipes to meet future heat transfer 
requirements covering a broad spectrum of transient operating conditions. 
These requirements are assumed to include high temperatures, high thermal 
flux, and reliable operation over long distances and small temperature drops. 
The model is intended to calculate the temperatures, velocities, and pressures 
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required to assess heat pipe transients. Transient operation characteristics of 
interest include the capability of a heat pipe to reach steady state operation 
following a change in operating conditions; the allowable magnitudes and rates 
of changes; the time required to respond to transient conditions; and the 
temperatures, velocities, and pressures encountered during transient operation. 
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Chapter II 

THEORY OF THE HEAT PIPE 


“[ The invention was to~\ cause absorption of heat, or in other words, the 
evaporation of the liquid to a point above the place where the condensation 
or the giving off of heat takes place without expending upon the liquid any 
additional work to lift the liquid to an elevation above the point at which 
condensation takes place.'' 


R. S. Gaugler 1944 [/2] 


2.1 HISTORICAL PERSPECTIVE 

The predecessor of the heat pipe was developed in the mid 1 9th century by 
Angier Perkins [6]. The device, known as the Perkins Pipe, consisted of an 
unwicked sealed container in which heat transfer was accomplished through 
evaporation and vapor transport to a cool condensing point. Circulation of 
liquid was accomplished through gravitational forces. This device was a 
wickless heat pipe, known in current terminology as the thermosyphon. 

The novel idea of using capillary forces to replace or supplement 
gravitational forces in circulating the working fluid is attributed to R. Gaugler 
in 1944 [12] . The device w r as a true heat pipe consisting of a sealed container 
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with capillary wick and working fluid. The invention was intended for use in 
the refrigeration industry but was not utilized. 

The development of space power systems some twenty years after Gaugler's 
patent introduced a viable application for the capillary heat transfer device. 
The term heat pipe was introduced by Grover and his co-workers at Los 
Alamos Scientific Laboratory in 1963[13] . Limited experimental work by 
Grover in 1964 using water as the working fluid followed by liquid sodium 
showed the great promise of heat pipes in serving as highly efficient conductors 
[1]. The theoretical basis for heat pipe operation and analysis was first 
formulated by Cotter in 1965 [2]. 


2.2 PHYSICAL OPERATING PRINCIPLES 

The term heat pipe refers to a large class of heat transport devices. The 
operating principles provide great flexibility in tailoring the heat pipe design 
to a specific application. The fluid dynamics and heat transfer principles are 
integral to the specific design configuration such that each unique design 
requires a dedicated analysis. However, the general operating principles can 
be described in terms of the simple representative configuration shown in 
Figure 3. 

A heat pipe generally consists of a long sealed container enclosing a working 
fluid and a capillary wick structure. The liquid and vapor have distinct flow 
paths, but share a common interface which allows them to communicate in all 
parts of the system. The quantity of working fluid is sufficient to at least 
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saturate the wicking material. The wick can be constructed of essentially any 
porous material: screen, gauze, sintered material, packed spheres, grooves in 
the container wall, perforated shields, etc. The heat pipe can use any working 
fluid that has a liquid and vapor phase at the operating temperature, that wets 
the wicking material, and is chemically compatible with the container and wick 
material. Proper selection of working fluid, container, and wick allow the heat 
pipe to be designed for specific temperature ranges: the cryogenic temperature 
range (0 to 150 K), the low temperature range (150 to 750 K), or the high 
temperature range (750 K and above). 
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A heat pipe generally has three distinct zones: the evaporator section, the 
adiabatic section, and the condenser section. The adiabatic section is optional 
and serves only to span the distance between the evaporator and condenser. 

There are three basic types of heat transfer/fluid dynamics problems to 
consider: 

1 . Vapor flow in the core 

2. Liquid flow in the wick 

3. Phase transition between liquid and vapor 

These three basic problems are coupled not only with each other but also with 
the specifics of the particular heat pipe design. Under static conditions with 
no heat flux, the liquid and vapor are under saturation equilibrium conditions. 
Thermal energy supplied to the evaporator is conducted across the container 
wall. Conduction and convection transport the heat to the liquid-vapor 
interface. The heat flux increases the liquid temperature which raises the 
liquid saturation vapor pressure. A pressure imbalance results between the 
liquid saturation vapor pressure and the pressure of the bulk vapor. If the 
pressure drop from the liquid saturation vapor pressure to the vapor phase 
pressure is larger than the pressure drop required for phase change, net 
evaporation will occur. The evaporating mass flux carries the thermal energy 
transferred to the liquid-vapor interface through the latent heat of 
vaporization. Mass injection and increased temperature raise the evaporator 
vapor phase pressure and establish pressure gradients that drive the vapor 
through the adiabatic region to the condenser. The pressure of the vapor in 
the condenser becomes higher than the liquid saturation vapor pressure. The 
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pressure imbalance results in net condensation on the cooler liquid if again the 
pressure drop is greater than that required for phase change. Since thermal 
energy is transferred in the form of latent heat of vaporization rather than 
sensible heat, high heat transfer rates can be achieved with low mass flow. If 
heat is transferred by high density, low velocity vapor the heat transfer is 
nearly isothermal due to the low vapor pressure gradient. 

The condensate is returned to the evaporator through capillary pumping in 
the wick. As liquid evaporates, the liquid-vapor interface recedes into the 
wick. Pores in the wick form capillary structures in which a curved meniscus 
forms as shown in Figure 4 due to the pressure difference between the liquid 
and vapor across the surface and the action of surface tension forces. The 
Laplace-Young Equation relates surface curvature to the pressure jump across 
an equilibrium liquid- vapor interfaced 14] 



where P v is the vapor pressure, P l is the liquid pressure, o is the surface tension 
coefficient, and R x and R 2 are the principal radii of curvature of the meniscus 
surface. The liquid contact angle is taken in Equation 1 as zero for simplicity. 
The effects of phase transition can be included to take account of 
non-equilibrium systems 

p v ~ p l+ AP pc = (2) 

where A P pc is the pressure difference due to the phase change transition. 
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As liquid recedes further into the wick, the radii of curvature decrease 
resulting in an increase in the liquid to vapor pressure jump and corresponding 
decrease in liquid pressure. The action in the condenser is the opposite process 
to that in the evaporator. Condensate tends to flood the wick such that the 
principal radii of curvature increase and the pressure difference between liquid 
and vapor decreases. If the wick is completely saturated such that the liquid 
forms a smooth flat surface across the wick at the liquid-vapor interface, then 
one principal radius of curvature is essentially infinite in length while the other 
is set by the dimension of the vapor space 

P,-P, + AP pc = -^- (3) 

where R v is the radius of the vapor space. A P pc is generally small relative to 
other pressure drops in the system and can be neglected. The quotient ajR v is 
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small relative to P v such that the pressure jump P v - Pi for the flooded wick 
can be neglected relative to other pressure drops in the system. The liquid and 
vapor pressures for a flooded wick can then be taken as being essentially equal 
as in the case of a completely planar interface. 

The variation along the length of the pipe of the radii of curvature coupled 
with the vapor pressure gradient results in a liquid pressure gradient lite'drives 
the liquid from the condenser to the evaporator. Typical pressure gradients in 
the liquid and vapor arc shown in Figure 5. These pressure gradients show 
that continuous circulation of the working fluid and corresponding heat 
transfer can be maintained. However, there are certain operating limits 
governed by fluid dynamics and heat transfer theory. Knowledge of these 
limits imposed by fundamental theory is required to assure heat pipe 
operation. 
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Figure 5: Liquid and Vapor Pressure ‘Gradients in a Heat Pipe 

a ) Negligible Gravity 
b ) Gravity Opposed 
Adapted from Ivanovskii et al. [4] 
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2.3 OPERATING LIMITS 


A heat pipe behaves like a structure of very high thermal conductance, but 
is constrained by the principles of fluid dynamics and heat transfer. These 
operating constraints depend additionally on thermophysical properties of the 
working fluid, container and wick, design configuration, and the coupling of 
the heat pipe to the environment for heat addition and removal. The 
commonly accepted limits discussed in the literature and shown in Figure 6 
include [4][15] : 

1. Viscous Limit 

2. Sonic Limit 

3. Entrainment Limit 

4. Capillary Limjt 

5. Boiling Limit 

As shown in Figure 6, the limits have overlapping regions of influence such 
that knowledge of all limits is required to assess a heat pipe for a specific 
application. 


2.3.1 Viscous Limit 

The viscous limit refers to vapor pressure losses due to friction forces. The 
viscous limit is generally associated with operation at low vapor pressures and 
temperatures when vapor density is low. The viscous limit is particularly noted 
in heat pipes with long condenser sections. Under these circumstances, the 
vapor pressure at the end of the condenser may be low enough to be effectively 
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zero in considering the maximum possible vapor pressure drop available to 
drive the vapor flow. The low possible pressure drop results in a low possible 
heat transport capability. The viscous limit is not necessarily a failure of the 
pipe to operate, but a limit on the heat transfer capability under specific 
conditions. 
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2.3.2 Sonic Limit 


The sonic limit refers to vapor flow attaining the speed of sound in some 
part of the pipe. This condition is generally a concern during startup when the 
vapor has low pressure and large specific volume. The vapor velocity at the 
evaporator exit is high even for low heat transfer rates and can reach sonic 
velocity. Mass flow rate is then choked at the point where sonic velocity is 
attained such that downstream pressure changes are not transmitted upstream. 
Since mass flow is choked, the axial heat transfer rate cannot be increased by 
improving heat removal in the condenser. In fact, improved condenser heat 
removal would actually prolong the sonic condition since the pipe would be 
prevented from heating up through storage of sensible heat. Large axial 
temperature gradients develop as the evaporator heats up faster than the 
condenser. The heat transfer associated with the sonic limit can be 
approximated as a function of vapor properties at the point of sonic velocity 
by [5] 


2 . 


son 


P v fyg A v ^ 

V2(l +y) 


(4) 


where p v is vapor density, c is the speed of sound corresponding to stagnation 
temperature, h fg is the latent heat of vaporization, A v is the vapor flow 
cross-sectional area, and y is the ratio of principal specific heats. 

The sonic limit ordinarily does not result in heat pipe failure but is a 
transition condition from low vapor pressure to normal operating pressure. As 
shown in Figure 6, sonic limit can prevent the heat pipe from reaching more 
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severe limitations that can disrupt heat pipe operation such as the capillary 
limit. 


2.3.3 Entrainment Limit 

The entrainment limitation is brought about by the capture of condensate 
in the wick by the counterflowing vapor. The process is similar to a body of 
water agitated by high velocity winds into waves which propagate until liquid 
is torn from the wave crests. In a similar manner, high velocity vapor in the 
heat pipe can tear liquid from the wick. The entrained liquid is carried with 
the vapor back into the condenser. Entrainment acts to increase the working 
fluid circulation rate. However, an increase in heat transfer occurs only if the 
liquid absorbs sensible energy prior to entrainment since phase change latent 
heat is not involved. Since sensible energy transport is much less efficient than 
transport of latent heat, the increased fluid circulation rate does not 
sufficiently increase heat transfer to meet the heat transport requirements. 
Fluid circulation increases until capillary pressure cannot accommodate the * 
flow rate. Formation of hot spots and dryout of the evaporator region results. 

The occurrence of entrainment is a result of vapor inertia forces overcoming 
liquid surface tension forces. Onset of entrainment is generally assumed to 
occur when the ratio of viscous forces to liquid surface tension forces, called 
the Weber number W e , is on the order of unity, where [4] 


We = 


2 i 

P v “ v * 
2nc 


( 5 ) 
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where u v is the vapor velocity axial component and X is the wavelength of a 
capillary wave. The capillary wave is a difficult to determine function of the 
liquid-vapor interface surface and the geometry of the surface pores. The 
entrainment limit heat transfer can be approximated as [5] 

Qem = (6) 

where r hp is the hydraulic radius of the wick pores at the liquid-vapor interface 
surface. 


2.3.4 Capillary Limit 

Continuous circulation of working fluid is required for steady state heat 
pipe operation. The liquid circulation is maintained by capillary pressure 
forces developed in the wick structure at the liquid-vapor interface. The 
capillary pressure forces A P cap must be sufficient to balance pressure drops in 
the vapor flow, the liquid flow, and phase transition 

AP cap >AP v + AP 1+ AP pc (7) 

As heat transferred by the pipe increases, the working fluid circulation rate 
also increases resulting in larger pressure drops in the flow paths. The 
maximum heat transfer capability is determined by the maximum capillary 
pressure force that can be attained by the wick. If the heat transfer is 
increased further, liquid flow rate cannot be increased to meet the higher 
demand. Evaporation will cause the liquid to recede into the wick resulting in 
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dryout and interruption of the working fluid circulation. The pipe will develop 
hot spots where liquid has dried out. The wicking limit is usually encountered 
at high vapor pressure when almost all pressure drop is in the liquid [15]. 

A simple expression can be derived relating the capillar}' heat transfer limit 
to operating conditions [5] 


AP 

Ax 


cap 


= “ F lQi 


cap 


( 8 ) 


where F ) is the liquid flow coefficient of friction, 0^ is the capillary heat 
AP 

transfer limit, and — — is the pressure gradient at the capillary limit. 

AX cap 

By assuming liquid and vapor pressures are equal at the end of the 
condenser, Equation 1 can be evaluated at the two extremes of the pipe 


- [/’y- ■P/] =[/>,•-/>,] =^- (9) 

x=0 x=L x=0 R 


where the assumption has been made that the capillary radii of curvature are 
equal, R { = R 2 = R. 

The capillary pressure drop limit over the length of the pipe occurs when the 
minimum capillary radii of curvature along the pipe are equal to the pore 
radius, R = R p , so that at the capillary limit 

AP a,p=jr (| 0 ) 


and Equation 8 can be written 


Qcap 
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where L is the pipe length and the frictional coefficient is given by 


F = H 

1 KA[hf g pi 


where A, is the liquid flow cross-sectional area, p, h fg , and p are the kinematic 
viscosity, latent heat of vaporization, and density, respectively, and K is the 
wick permeability given by 


K = 


2er[ 

fl Re l 


with £ the flow direction porosity, r h the hydraulic radius, and f Re, the friction 
factor- Reynolds number product for laminar duct flow. 


2.3.5 Boiling Limit 

The boiling limit refers to the occurrence of nucleate boiling in the wick 
structure. While many convection systems are enhanced by nucleate boiling, 
formation of vapor bubbles can block flow passages in wicks and cause 
localized hot spots. The resulting increased circulation rate can exceed the 
capillary capabilities of the pipe and disrupt operation. Boiling limit is highly 
dependent on the working fluid thermophysical properties. A working fluid 
with low thermal conductivity such as water may easily boil in the wick 
structure, while a high conductivity working fluid such as liquid metal may 
efficiently conduct energy to the liquid- vapor interface and resist boiling[16] . 
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The extensive variables involved with boiling processes make the boiling limit 
the most difficult to accurately predict. 
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Chapter III 

HEAT PIPE DYNAMICS 


“ The quantitative details of the startup process are of minor interest." 

T. P. Cotter 1965 [2] 

Heat pipe dynamics are difficult to predict due to the complex fluid 
dynamics and heat transfer phenomena and the wide range of variables 
involved. These variables include pipe configuration, working fluid fill, initial 
temperature, the rate of change of heat addition to the evaporator, and 
condenser heat ‘sink conditions. The important dynamic processes are also a 
function of the heat pipe operating state at each point during the transient. 
Transient operation can be considered in terms of three general regimes: 
startup; operating transients; and shutdown. The following sections discuss 
the physics of these regimes and the status of research in each area. 
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3.1 STARTUP 


Heat pipe startup refers to the process of bringing a heat pipe from an 
initial static condition to steady state operation. The heat pipe working fluid 
may initially be in a frozen or liquid state. Heat pipe startup concerns include 
the conditions and time necessary to bring a heat pipe to operation, the 
conditions that hinder startup, and the dynamic behavior during the startup 
period. 


3.1.1 Startup Physics 

There are various means to attain heat pipe startup from either liquid or 
frozen state. One reliable method is to isothermally warm or cool the entire 
pipe to the desired operating temperature and then gradually increase the 
input heat flux. Isothermal conditioning of the entire heat pipe may not be 
practical, particularly for the long heat pipes required for future applications 
such as for reactor core heat removal. Other methods of startup can be. 
grouped according to the three general startup modes as shown in Figure 7. 

Uniform startup corresponds to Figure 7(a). Uniform startup occurs when 
the density of the vapor is high at the initial temperature. The working fluid 
is completely melted. Heat applied to the evaporator results in vapor pressure 
gradients that initiate working fluid circulation. The high vapor density and 
pressure result in a high rate of heat transfer with a low vapor pressure drop. 
Since axial temperature drop is a function of vapor pressure drop, the axial 
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temperature gradient is also small. Heat pipe response is smooth and relatively 
quick. 
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Frontal startup is shown in Figure 1(b). Frontal startup occurs when the 
vapor density at the initial temperature is low such that vapor flow is free 
molecular. Free molecular flow refers to the molecular mean free path 
exceeding the diameter of the vapor channel. At low vapor pressure the 
molecular mean free path is large so that interaction between vapor molecules 
is rare. The nature of flow is determined mainly by collisions of vapor 
molecules with the pipe walls. Evaporator heating raises vapor pressure in the 
evaporator so that vapor begins to flow as a continuum when the molecular 
mean free path becomes small compared to the vapor channel diameter. 
Vapor in the condenser is still low density and free molecular, so that there is 
little axial heat flux. A transition region exists between evaporator and 
condenser. Vapor compressibility effects become important since vapor flow 
may reach supersonic velocities. Since axial temperature drop is determined 
by vapor pressure drop, evaporator temperatures will initially be much higher 
than in the condenser. Condensation of liquid drops occurs in the vapor space 
during the rapid expansion through ‘the transition region. The transition 
regime occurs at intermediate mean free paths such that collisions between 
vapor molecules and collisions with the wall are of the same order of 
magnitude such that both types of interactions determine the flow. The 
entrainment limit may be reached due to the viscous stress from high vapor 
velocity tearing liquid from the wick. 

The heat pipe temperature and vapor pressure increase since the sonic limit 
restricts the flow of latent heat and forces the retention of sensible heat. At 
high vapor pressure the molecular mean free path is small so that vapor 
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molecules interact frequently and determine the nature of flow as continuous. 
The condenser temperature will then increase until the pipe becomes almost 
isothermal. Startup is completed as with a uniform startup. 

Vapor flow regimes are determined by vapor temperature and pressure. 
For a sodium heat pipe, Ivanovskii et al. found free molecular flow from the 
melting point to 250 °C, transition flow from 250 to 400 °C, and continuous 
flow above 400 °C [4], Criteria were suggested for bounding the flow regimes 
based on the Knudsen number: 


^ c 

where L p is the molecular mean free path and L c is a characteristic dimension 
for vapor flow in a heat pipe such as the vapor channel diameter D v . The 
suggested boundaries are given as: 

Free Molecular Flow 
Transition Flow 

— — < 0.0 1 Continuum Flow 

D v 

Startup may fail to occur if initial vapor pressure is low and the thermal 
resistance between condenser and heat sink is low. Low thermal resistance 
maintains a good flow of heat out of the pipe so that the pipe does not retain 
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L p 

0.01 < — -< 1.0 
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the energy to raise the overall pipe temperature. Low vapor pressure results 
in compressibility effects and possible sonic conditions. Startup may fail due 
to the entrainment limit. The pipe may be successfully started from low vapor 
pressure if there is a high thermal resistance at the sink-container interface. 
While vapor flow may reach the sonic limitations, the condenser temperature 
will gradually rise and raise the vapor pressure so that normal operation is 
achieved. 

The working fluid may originally be in a frozen state. Heat introduced to 
the evaporator progressively melts the working fluid to. the empty vapor space. 
Evaporation occurs to essentially a vacuum such that free molecular flow 
occurs above the melting and solid working fluid. Vapor molecules that 
condense on the frozen fluid will then refreeze. Heat pipe startup can be 
prevented if evaporation and subsequent condensing and freezing of working 
fluid occurs at a faster rate than the melting for replenishing the supply. The 
rate of heat input is an important consideration in achieving startup, as well 
as maintaining the condenser insulated so that heat is kept in the pipe to melt 
the working fluid. 

A heat pipe pressurized with a noncondensible gas may also undergo frontal 
startup as shown in Figure 7(c). A gas loaded heat pipe is useful as a 
temperature regulator or as a thermal diode. The presence of noncondensible 
gas is also known to aid the startup of a heat pipe from a low initial vapor 
pressure. As working fluid evaporates, a diffuse vapor-gas transition region 
occurs creating a thermal front and forming a compressed gas plug in the 
condenser.[l 7] . Additional evaporation further compresses the gas such that 
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eventually only part of the condenser vapor space is plugged by gas. The 
length of this inactive condenser zone is a function of the vapor pressure 
dependence on temperature. Since the vapor is in a near saturation condition, 
the vapor pressure is directly related to the vapor temperature. The axial 
temperature distribution therefore provides an accurate description of the 
corresponding vapor conditions[18] . A qualitative axial distribution for 
frontal startup of a heat pipe with noncondensible gas is shown in Figure S. 



Figure 8: Frontal Startup Temperature History of a Gas-Loaded Heat 

Pipe 

Adapted from Bystrov and Goncharov [18] 
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The heat pipe is initially at a uniform temperature represented by curve 1. 
A uniform applied heat flux causes a uniform warming of the evaporator as 
represented by curve 2. The noncondensible gas overpressure is assumed to 
suppress evaporation during the initial warming period. As the temperature 
in the evaporator reaches curve 3. intense evaporation begins. The mass 
addition sweeps the noncondensible gas into the condenser. Sonic velocity is 
reached at the evaporator-condenser boundary. The vapor velocity 
corresponding to curve 3 is slowed to subsonic after the evaporator exit due to 
the buffering effect of the noncondensible gas. The velocities corresponding to 
curves 4 and 5 reach supersonic conditions at the evaporator exit as shown by 
the sharp temperature drop. Subsequent pressure and temperature recovery 
occur at the noncondensible gas front. 

As vapor pressure increases the noncondensible gas front is compressed 
further into the condenser. Supersonic flow occurs over the length of the active 
condenser as shown by curves 6 and 7. Non-isothermal conditions during this 
period of sonic limited operation can be very pronounced. The shock front 
adjustment from supersonic to subsonic velocity gradually moves to the end 
of the condenser. The shock front reflects back to the evaporator as 
represented by curves 8 through 10. The vapor speeds during this process are 
sonic at the choked evaporator outlet, supersonic in the condenser up to the 
shock, and subsonic after the shock. As shown by curve 11, continued 
warming raises the vapor pressure and eliminates the shock conditions. While 
initially a large axial temperature drop persists, the heat pipe eventually 
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reaches near isothermal conditions as shown by curves 12 and 13. High heat 
transfer rates occur during this period of near isothermal conditions. 


3.1.2 Experimental Startup Studies 

Experimental studies were performed by Ivanovskii et al. on the 
temperature distribution along a heat pipe during startup [4], The 
experiments were performed with sodium heat pipes with and without 
noncondensible gas added to the vapor space. The entire inventory of working 
fluid was solid at the beginning of the experiments. The experiments studied 
the frontal startup mode. Temperature measurements made in the vapor core 
are shown in Figure 9. 

Curves 1 through 5 show that axial heat transfer to the condensation zone 
is negligible during initial heating of the evaporator. As the temperature 
increases, the flow of vapor molecules from the evaporator to condenser 
increases. The authors state that curves 5 and 6 represent free molecular and 
transition flow regimes. The flow regime can be seen to be a function of local 
axial conditions. 

Curves 6 through 9 represent operation at the sonic limit with supersonic 
flow in the condenser due to low vapor pressure and high condensation rate. 
A shock front occurs corresponding to the abrupt temperature drop at the 
transition from supersonic velocity to zero velocity at the end of the condensing 
zone. Non-isothermal conditions during sonic limited operation can be seen 
to exceed 300 °C. Continuum vapor flow begins with curve 8. 
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Curves 10 through 14 represent steady operation below the sonic limit. The 
authors point out that the degree of isothermal operation depends on the vapor 
pressure. The extent of isothermal operation increases with decreasing ratio 
of heat transferred to the sonic limit. Heat pipe operation is essentially 
isothermal for the large vapor pressures represented by curves 12 through 14. 
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The time dependence of the startup was not reported, although the authors 
state without definition or elaboration that the warmup time was ninety 
minutes. 

Experiments were performed at Los Alamos Scientific Laboratory on a heat 
pipe with molybdenum container using lithium as the working fluid[19][20] . 
The intent of the program was to address questions pertaining to the operation 
and performance characterization of high length-to-diameter ratio heat pipes 
under conditions of high radiation loading. Startup times and power handling 
capability, shutdown behavior with radiation sink temperatures below the 
freezing point of the working fluid, and accuracy of presently used 
performance models were specific questions addressed. 

Starting from a uniform temperature of approximately 300 K, the heat pipe 
evaporator was brought to a uniform temperature of 1350 K over a period of 
forty minutes. The temperature was then held constant by adjusting the input 
power as the melt front progressed along the heat pipe. Total time to bring the 
heat pipe to temperature was approximately three hours. Peak power 
throughput was 15 kW. The temperature history of a typical condenser 
thermocouple (located 106 cm from the evaporator end) is shown in Figure 10. 
The interesting characteristic is the sharp thermal front between the melted 
and solid regions of the pipe. 

Transient operation of low temperature heat pipes has been studied 
experimentally by Colwell et al. for several years[2l3E223[233 . Experimental 
studies concerned a heat pipe with circumferential and central slab wick using 
Freon 1 1 as the working fluid. Testing was performed in temperature ranges 
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Figure 10: Startup Temperature History of Experimental Heat Pipe 

Adapted from Merrigan et al. [19] 


from ambient room conditions to well above the working fluid critical point. 
The heat pipe was found to operate very smoothly at modest heat transfer 
rates with temperatures well below critical. Abrupt changes in coolant flows, 
heating rates, and inlet coolant temperatures were easily accommodated. 
Liquid temperatures were nearly uniform in the slab wick and only small 
gradients existed in the vapor. Measured vapor pressures were very close to 
saturation pressures corresponding to measured temperatures. Steady state 
was generally reached within ninety minutes for the specific configuration and 
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testing performed. Testing conducted at temperatures near the critical 
temperature identified more dynamic conditions in the heat pipe than at the 
lower temperatures. Partial and complete drying of the wick were 
encountered. Experiments were also performed to study the rewetting 
capabilities of the wick after dryout. Under some conditions steady operation 
was achieved even with a portion of the evaporator wick dried. Rewetting 
could sometimes be accomplished depending on the change in boundary 
conditions. Rewetting was found to be very difficult to accomplish after the 
wick completely dried. A significant lag time was encountered in rewetting the 
wick. 

Mathematical models were formulated to describe the experimental 
results[24] . The model for a fully wetted wick was based on the assumptions 
that heat is transferred by conduction through the container wall and liquid 
saturated wick, thermal resistances associated with phase change are 
negligible, and vapor space temperature is uniform. A finite difference code 
using the Alternating Direction Implicit method was used to solve the 
conduction equation. The vapor temperature was determined from a simple 
heat balance at the liquid-vapor interface. Under saturated wick conditions 
the transient conduction equation with variable properties predicted 
performance accurately. A lumped parameter one dimensional momentum 
equation was used for cases in which the wick was not completely saturated. 
The mathematical model did not agree with experiments for transients 
involving a partially dried wick. General trends and end steady state 
conditions were qualitatively predicted, but with large temporal errors. A 
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calculation time step of five seconds was used, although the authors state that 
a time step of up to twenty seconds could be used. 


3.1.3 Analytical Startup Studies 

The first apparent study of heat pipe transients was reported by Cotter in 
1967[25] . A one dimensional analysis was performed of a heat pipe at 
uniform initial temperature 7, which is heated to steady state operating 
conditions at temperature T 2 and heat flux Q 2 . A uniform heat flux O in is 
supplied to the evaporator during startup. Q m may vary with time depending 
on the type of startup to be considered. A slow startup consists of increasing 
Qtn t0 Q2 so slowly that the pipe essentially passes through a series of steady 
states. A rapid or pulsed startup occurs for Q in = 0 at time t = 0 and 
Qm = Qi f° r firne t > 0. By considering the pipe as a one dimensional object 
with temperature T a function of axial location x and time t. Cotter derived 
energy equations for the evaporator and condenser regions 



+ C-^— — ^‘ n 
dt ~ l e 

_ -Qi 
lc 


T- T x 


for the evaporator 
for the condenser 


(13) 


where the effective thermal conductivity in the axial direction is k eff , A is the 
pipe cross-sectional area, C is the heat capacity per unit length, and l c is the 
condensation zone length. 
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Equation 13 suggests two different characteristic times. The first, 
CIJ,T 2 - T x )IQ 2 , is related to the heat input required to raise the pipe 
temperature from initial value T ] to final value T 2 . This first characteristic 
time is given as on the order of 10 to 100 seconds. 

The second characteristic time, CQIk ef/ A, is related to the axial heat 
transport rate. This time is inversely proportional to the effective thermal 
conductivity. At low vapor temperature and corresponding low vapor 
pressure, the molecular mean free path is large compared to the pipe diameter 
so that heat transfer is primarily due to conduction in the pipe wall and liquid. 
The second characteristic time under these conditions is given as on the order 
of 10 3 seconds. As the vapor pressure increases, convective heat transfer of 
latent heat becomes of comparable magnitude to thermal conduction. As the 
vapor temperature reaches a temperature range of 100 to 200 °C, the transition 
to continuum vapor flow is completed. The effective thermal conductivity 
increases rapidly such that the second characteristic time is typically less than 
one second. 

Cotter performed a first approximation analysis based on the relative 
magnitudes of the characteristic times with the results shown in Figure 1 1 . 
Let T c be the temperature at which continuum vapor flow is established. Due 
to the large characteristic time, heat transfer is neglected wherever the 
temperature is below T c . Wherever the temperature in a region exceeds T c , 
heat transfer is so rapid that the entire region has uniform temperature. By 
considering the bounds set by rapid startup and slow startup Cotter solved a 
system of equations for startup from various initial temperatures T v Startup 


38 


from a high vapor pressure in which 7, > T c results in a uniform startup and 
the hot zone extends over the entire length of the pipe. The uniform startup 
is plotted as line segment AC in Figure 11. Solution for the rapid startup is 
plotted as line BC in Figure 1 1. Any startup program lying within the region 
ABC is expected to be successful. 



Startup from low vapor pressure with 7, < T c will have no initial axial heat 
transport as represented by point D in Figure 11. The evaporator heats up 
until 7 = T c at point E. The sonic limit line of curve 4 is followed above T c to 
point F. A slow startup will follow curve FC until steady state operation is 
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reached at point C. A rapid startup follows a slightly peculiar path. The 
heated zone reaches the end of the condenser on the interval GH , where G is 
the intersection of the sonic limit curve and curve GC. At point H, O son = 0 2 , 
but T 2 has not been reached. The heat pipe continues to heat up such that the 
path H—I—C is followed. The startup regime for a pipe at initial 
temperature below T c is then the curve DEF and the region FGH1C. 

The heat transfer variation with time was expressed by Cotter using a one 
dimensional approximation. Heat supplied to the evaporator was expressed 
as the exponential function 

Qin = Q2I 1 “ ex P( “ '/ 'o)] ( 1 4 ) 

where t 0 = Cjk efJ -. With t Q = CIXT 2 - T x )jQ 2 , the variation of heat transfer with 
time is given by 



The qualitative behavior of Equation 15 for a range of heating rates is shown 
in Figure 12. 

Researchers have followed Cotter's startup analysis with approaches that 
have ranged from the simple to the complex. The simple end of the spectrum 
includes the adaptation of generalized thermal analyzer programs to include 
heat pipe transients[26][27] . The intermediate range of complexity includes 
lumped parameter approaches, while the high end of complexity involves 
solution of the full system of governing equations. 
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Figure 12: Analytical Heat Pipe Startup as a Function of Time 

( 1 ) Rapid Startup 
( 2 ) Intermediate Startup 
( 3 ) Slow Startup 
Adapted from Cotter [25] 


Published lumped parameter transient studies include two studies of 
thermosyphons. A thermosyphon is similar to a heat pipe in that heat is 
transferred in a closed container through the latent heat of vaporization of the 
working fluid. The two devices are dissimilar in the mechanism for returning 
the liquid phase to the evaporating region of the device. A thermosyphon uses 
gravitational force to accomplish liquid return while a heat pipe uses capillary 
forces that are either opposed by gravity, supplemented by gravity, or are 
gravity free. The similarities between devices make analysis of thermosyphons 
of interest in the study of heat pipes. 


41 




A simple lumped parameter formulation for thermosyphon transients was 
developed by Dikiy et al.[28] . The thermosyphon considered was actually 
simplified from the usual configuration since the vapor flow paths and liquid 
flow paths were assumed to be separate and isolated so that communication 
between the countercurrent liquid and vapor was not considered. Heat and 
mass balances were formulated for liquid, vapor, evaporator wall, and 
condenser wall. Empirical relationships were used in defining heat transfer 
coefficients for the condensation and evaporation processes. The result of 
derivations was a first order nonlinear differential equatio'n for temperature 
as a function of time. Results of the numerical integration of the governing 
differential equation show transient processes that occur smoothly over a 
period of approximately thirty seconds. 

Reed and Tien developed a lumped parameter model of a thermosyphon 
that included the effects of liquid-vapor interaction[29] . Control volume 
equations were developed for conservation of mass, momentum, and energy in 
the liquid; mass, momentum and energy in the vapor core; and overall system 
mass and energy. Auxiliary empirical and analytical relations were used for 
specifying wall and interfacial shear stresses and the liquid film heat transfer 
coefficient. The lumped parameter model compared well with steady state 
experimental data. Transient experimental data were not available for 
comparison. The primary significance of the approach was that a single model 
was found to be able to predict most operating parameters and operating 
limits, making it useful for inclusion in a systems analysis formulation. 
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A lumped parameter approach was used by Beam et al. to investigate heat 
pipe response to pulsed startups[30][31] . A pulsed startup refers to applying 
a step increase in heat flux to a static heat pipe. The authors address the 
questions of whether a heat pipe will reach normal operation after applying the 
pulse, and the time required to rewet a wick that partially drys during the 
startup. The time variation of heat pipe temperature was expressed as an 
exponential function. The authors also developed a model to describe the 
dryout and rewetting of wicks. Two types of pulsed startups were considered. 
In the first the mean temperature of the transport section was maintained 
constant at the initial temperature so that thermal storage in the system was 
negligible. This test addressed the ability of the heat pipe to respond 
immediately to pulsed heating with no thermal storage. In the second, the heat 
pipe temperature was allowed to increase so that thermal storage was 
appreciable. The lumped parameter model was compared to experimental 
data from a copper-water heat pipe with a homogeneous screen wrap wick. 
Three experimental responses to pulsed startup were identified. The wick in 
the evaporator may accommodate the pulse heat flux without drying; partial 
dryout may occur followed by rewetting; or dryout may occur without 
rewetting. The authors conclude liquid in a screen wick can respond almost 
instantaneously to a change in pressure gradient. The heat pipe was found to 
be able to transport a pulsed heat load equal to the capillary pumping limit 
without drying and without thermal storage. 

Experimental and analytical work addressing startup with and without 
noncondensible gas was performed by Bystrov and Goncharov [18]. Two 
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analytical approaches were attempted. In the first, conservation of mass, 
momentum and energy equations for vapor were written as one dimensional 
averaged over the vapor channel cross-section. Phase change from evaporation 
and condensation was based on kinetic theory. The vapor phase equations 
were solved separately for the evaporator and condenser zones with 
appropriate matching conditions at the boundary between zones. The liquid 
phase solution method depended on whether the working fluid was initially 
frozen or solid. The non-steady heat conduction equation was used for the 
solid case with a heat balance to describe melting. The one dimensional 
equations of motion for a flux of variable mass were used for the liquid phase. 
The area of liquid flow was taken as constant. Heat supply and rejection were 
through radiation boundary conditions. 

The resulting numerical solution routine required extremely small time 
steps on the order of 10~ 5 to 10~ 7 second for stability. The equations governing 
the liquid phase had better stability characteristics such that time steps of 
10~ 2 second could be used. The authors note that these characteristics have 
advantages when the vapor flow is subsonic. The vapor can be taken as 
uniform due to large time scale differences. 

The small time steps required for vapor were prohibitive for a complete 
transient analysis. The authors developed a revised approach consisting of 
lumped parameter heat balances that avoided calculation of detailed 
temperature distributions. The kinetic model of interphase mass transfer was 
also eliminated to avoid the stability problems. The authors solved the heat 
balance equations, one dimensional conservation of mass, momentum and 
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energy for the vapor, conservation of mass and momentum for the liquid, and 
supplemental equations to maintain operation within known operating limits. 

The comparison of analytical and experimental results for a heat pipe with 
noncondensible gas is shown in Figure 13. The startup exhibits frontal 
behavior as noncondensible gas is swept into the condenser. The region 
furthest in the condenser has a large lag time showing the low diffusion of 
vapor into the gas. Initially the evaporator warms uniformly while there is 
little axial heat transfer since the noncondensible gas suppresses evaporation. 
A substantial axial temperature gradient in excess of 400 °C occurs due to 
sonic conditions at the evaporator outlet. 

A detailed computer model of transient heat pipe performance has been 
under development at Los Alamos Scientific Laboratory[32] . The model is 
intended to include startup from frozen state with provision for free molecule 
flow through steady state operating conditions. Vapor flow is modeled as one 
dimensional and compressible with friction and mass addition from the wick. 
The Kachina method is used for solving the governing equations. The authors 
state mass addition from the wick is the most difficult aspect to simulate. 
Mass transfer rate based on kinetic theory was found to produce high effective 
mass transfer coefficients that required very small time steps for explicit 
integration. This problem was circumvented by using a simplified wick model 
for computing temperature and saturation pressure at the wick surface. Mass 
transfer rates based on this locally computed wick surface temperature are 
used in the vapor core solution and later in the full wick solution. The wick 
model (liquid flow) is solved using the SIMPLER method. Freezing and 
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Figure 13: Comparison of Analytical and Experimental Startup Data 

a ) Experimental Temperature Fields 
b ) Experimental and Analytical Data 

1 . Evaporator Surface Temperature 

2. Condenser Surface Temperature 

3. Condenser Average Temperature 

4. Evaporator Exit Temperature 

5. Heat Transfer 

Adapted from Bystrov and Goncharov [18] 


thawing are treated by an enthalpy balance. Radial temperature of the 
wick/groove/wall structure is solved implicitly at each axial location. Axial 
direction integration is performed explicitly. Implicit radial integration is used 
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due to the high radial conductances and the small radial dimensions. The 
model was run for ten computational steps on a Compaq Deskpro II 
computer. Execution on a Cray computer was expected to require thirty 
minutes of processor time per hour of actual time simulated. 

Analytical work by Peery and Best was performed to address rapid 
transients in a heat pipe[33] . The heat pipe simulated was a 200 cm long 
rectangular aluminum container with a homogeneous slab screen wick and 
water working fluid. Continuity and momentum were solved as one 
dimensional while energy was solved as two dimensional. The computations 
were fully implicit in time with property variations a function of temperature. 
The liquid and vapor flows were solved individually with coupling through a 
kinetic phase change model. The authors began computations with a time step 
of 1 x 10 -4 seconds and increased the time step to 6.4 x 10 -3 seconds. The 
model broke down for time steps larger than 6.4 x 10 -3 seconds. At three 
seconds into the transient, the condensation rate increased to become 
approximately equal to the evaporation rate. Convergence then became a 
problem and the transient was terminated. The code was expensive to use 
since one second of transient required one hour of CPU time on a CDC Cyber 
825 computer. The authors conclude that the kinetic theory phase change 
model caused the requirement for small time steps and resulting computational 
expense. 
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3.2 OPERATING TRANSIENTS 

Operating transients involve power changes during nominal operating 
conditions. There are several transient conditions of interest. A heat pipe may 
operate at low power standby mode but be required to quickly respond to 
reach full power operation. A heat pipe may also be required to perform load 
following functions depending on changing conditions at the power conversion 
system. Reactor control also may require changes of heat generation that 
would have to be handled by the heat pipe system. 

The physics of interest for these situations depends on the initial condition, 
the rate of change and magnitude of input power, and condenser conditions. 
A heat pipe operating in standby mode at low vapor pressure may encounter 
the sonic limit while powering up. Large temperature gradients and complex 
vapor dynamics would develop. After passing through sonic limited operation, 
or for a heat pipe initially operating at high vapor pressure, the capillary, 
entrainment, and boiling limits are of concern. 

While physical processes are similar to those encountered during a startup, 
very limited work has been performed to address issues or characterize 
transient performance. The previously discussed study performed by 
Ivanovskii et al. and shown in Figure 9 included some work in the operating 
transient regime. Curves 11 through 14 are the axial temperature profiles at 
different axial power levels. These curves show the temperature characteristics 
of a quasisteady process since the input heat flux was intentionally varied 
slowly with time. The issues concerning more dynamic processes remain to be 
resolved. Issues requiring resolution involve the allowable rate of change of 
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input heat flux, the magnitude of input heat flux, the time required for the heat 
pipe to respond to a transient, and the capability to buffer small scale 
transients that occur at the heat sink and source. 


3.3 SHUTDOWN 

The dynamics of heat pipe shutdown are important since subsequent restart 
depends on the final shutdown condition. Shutdown may be to a liquid or 
frozen state. Depending on the means of attaining shutdown, wick depriming 
may occur due to non-uniform conditions. Shrinkage voids may occur as the 
liquid begins to freeze. Shrinkage may also cause tearing of the delicate 
capillary wick structure. The tear would result in a larger capillary surface and 
corresponding reduction in capillary pumping capability. 

These issues are illustrated in the single identified study that addresses heat 
pipe shutdown dynamics [19] [20]. Shutdown tests were performed at Los 
Alamos Scientific Laboratory on a high temperature heat pipe with annular 
wick and lithium working fluid. A temperature profile taken 8 minutes 30 
seconds after power shutoff is shown in Figure 14. Lithium near the 
evaporator end of the condenser has reached the solidification point while 
lithium in the evaporator is almost 240 K above the melting point. An 
attempted restart of the pipe after complete lithium solidification failed due to 
voids in the annular flow region. The authors concluded the voids were a 
result of shrinkage of the working fluid during freezing due to non-uniform 
cooldown rates. Subsequent tests were performed in which a more uniform 


49 



cooldown rate was attempted. The temperature difference between the hottest 
and coolest point at onset of solidification was reduced to approximately 140 
K. Subsequent restart tests were performed successfully. The authors 
concluded the critical issues include the axial location where freezing first 
occurs during shutdown, and the ability of the wick structure to accommodate 
working fluid shrinkage in the region furthest from the condenser and liquid 
pool following establishment of the freeze plug. The restart capability of a 
given radiation coupled design was found to be strongly dependent on the 
external geometry and thermal environment of the heat pipe as well as on the 
heat pipe internal design itself. 
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Chapter IV 
THE SOLUTION 


“The simple derivation is left as an exercise for the interested reader 

Engineering Textbooks 


4.1 THE SCOPE 

As previously discussed, most research in heat pipe dynamics has focused 
on the startup transient. An experimental investigation has also been 
performed for the shutdown transient. Only modest efforts have been 
expended to date that consider operational transients. These efforts generally 
address two issues: the capability to achieve nominal operation under specified 
conditions, and the resulting axial temperature distribution. The processes are 
assumed to occur quasisteadily such that time dependent characteristics are 
not considered. Despite these research efforts, a comprehensive analytical 
model does not exist that captures the physical processes while being 
economical enough for required computations. 
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The objective of this study is to develop an analytical model of the liquid 
phase of a high temperature heat pipe. The model is intended to be coupled 
to a vapor phase model for the complete solution of the heat pipe problem. 
The mathematical equations are formulated consistent with physical processes 
while allowing a computationally efficient solution. The model simulates time 
dependent characteristics of concern to the liquid phase, including: input, 
phase change, and output heat fluxes, liquid temperatures, container 
temperatures, liquid velocities, and liquid pressures. 


4.2 THE CONFIGURATION STUDIED 

This research considers requirements for high temperature, high thermal 
flux operation with heat transport over long distances and small temperature 
drops. These requirements indicate the need for a heat pipe with a high 
performance composite wick. A composite wick consists of small pores at the 
liquid-vapor interface and large pores in the direction of bulk liquid flow. 
Recall that the maximum capillary pumping capability is limited by the largest 
capillary surface at the liquid-vapor interface as shown by Equation 1 . The 
small pores at the liquid-vapor interface establish a high capillary pumping 
capability that drives the circulation of working fluid. The larger pores in the 
direction of flow provide a low resistance flow path for the liquid. 

The heat pipe selected for this study uses an annular wick configuration as 
shown in Figure 15. The annular wick is an excellent configuration for liquid 
metals, is popular for experimental work, and involves geometry easily 
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modeled in cylindrical coordinates [15,19,20]- The wick configuration consists 
of several layers of fine pore screen pressed together and concentrically 
installed in the pipe to form an annular flow channel as shown in Figure 16. 
The annulus forms the low resistance flow channel while the fine pore screen 
tube forms the high pumping capability boundary between the iiquid and 
vapor regions. 



INPUT 

Figure 15: Heat Pipe with High Performance Composite Annular Wick 


The capillary wick introduces some analytical although not conceptual 
difficulty. A thin screen tube consisting of several layers of pressed screen 
forms a tortuous flow channel with no slip conditions on the surface of the 
screen wires. Explicit analysis of the flow including consideration of the large 
number of wire surfaces involved would be difficult if at all practical 
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considering the goal of performing transient analysis. A macroscopic approach 
such as Darcy's Law for flow in a porous medium could be used to alleviate 
the complexity. However the axial flow in the pressed screen is expected to be 
small relative to axial flow in the open annular channel. Indeed, the very 
purpose of the annular composite wick configuration is to provide a low 
resistance flow channel in the annulus for circulation of liquid. The screen 
tube is used only to provide a fine capillary surface at the liquid-vapor 
interface. The screen tube can then be approximated as a single layer of screen 
with a characteristic pore size. Flow normal to the screen is restricted by a 
characteristic specified porosity, while tangential flow is governed by 
alternating free slip and no slip conditions. The free slip condition provides the 
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tangential shear communication between the liquid and vapor phases. The 
wick is approximated then as an infinitesimally thin porous surface. 


4.3 THE LIQUID PHASE MATHEMATICAL FORMULATION 
4.3.1 Assumptions 

The governing equations to be presented below will be reduced by 
recognizing unique characteristics of the system to be analyzed. The system is 
assumed to operate in an essentially gravity free environment. Other body 
forces such as those due to translational and rotational acceleration will also 
be neglected. 

The heat pipe is assumed to be installed with azimuthal symmetry in the 
particular application. The thermal energy source adds heat uniformly around 
the evaporator circumference while the heat sink extracts heat uniformly 
around the condenser circumference. The resulting uniform conditions in the 
azimuthal direction imply there are no velocity components or gradients in the 
azimuthal direction. 

Liquid flow in a heat pipe wick is generally low speed, laminar, and 
incompressible. Since the liquid flow field is relatively simple and the 
cross-stream dimension of the liquid region is much smaller than the 
streamwise dimension, the cross-stream (i.e., radial) pressure gradient is 
assumed to be negligible. Also since the flow is low speed and incompressible, 
viscous energy dissipation can be neglected. 
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In summary, the assumptions to be used for reducing the governing 
equations are: 

1. Negligible body forces: f r =f x = 0 

2. Azimuthal symmetry 

a. No azimuthal flow: u- = 0 

b. No azimuthal gradients: = 0 

d(p 

3 p 

3. Negligible radial pressure gradient: — = 0 

dr 

4. Negligible viscous dissipation 

5. Incompressible liquid: -^-=0 

6. Laminar liquid flow 

Thermophysical properties are considered to be uniform functions evaluated 
at a reference temperature. The reference temperature is assumed to change 
slowly with time such that properties can be treated as constants in the 
governing equations. 

4 . 3.2 Governing Equations 

Temperatures of the heat pipe container and the temperatures, pressures, 
and velocities of the liquid phase are required to describe the dynamic behavior 
of the liquid phase of a high temperature heat pipe under transient operating 
conditions. The governing equations describing the behavior of these 
parameters consist of the conservation of mass, momentum, and energy. 
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4.3.2. 1 Conservation of Mass 

The full condition governing conservation of mass of an incompressible fluid 
is 

iiL + iL + J_i2L + _L =0 (16) 

ax dr '• d<f> r 

where u is the velocity in the axial (x) direction, v is the radial (r) velocity, and 
w is the azimuthal direction (<f>) velocity. Imposing the assumption of 
azimuthal symmetry reduces Equation 16 to the familiar two dimensional form 

_^L + 4l + _L = o (17) 

dx dr r 


4.3. 2.2 Conservation of Radial Momentum 

Full conservation of radial momentum for an incompressible fluid with 
constant properties is given by 


2 2 2 
dv dm dv . 1 dwv w . v _ 

dt dx dr r d<j> r r 

i dp ( a 2 v a 2 v i a 2 v i dv 2 dw 

P dr +V \ax 2 • dr 2 r 2 d(f> 2 r dr r 2 



(18] 


where t is the time, p is density, P is pressure, and v is kinematic viscosity. 
By assuming negligible radial pressure gradient, the conservation of radial 
momentum equation is replaced by 


-EZ_=0 (19) 

dr 
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Under this approximation the radial momentum equation is not used to define 
liquid dynamics so that pressure is determined from the continuity equation 
and the axial momentum equation. 


4.3.2.3 Conservation of Axial Momentum 

Conservation of axial momentum for an incompressible fluid with constant 
properties is given by 


du cn?_ 
dt dx 


duv _1_ duw _wv 

dr r d<f> — 


J_ d_P , f d^u d^u_ 1 d 2 u , 1 du \ 

\dx 2 5r 2 r 2 d<t> 2 r dr J 


+ fx 


( 20 ) 


Assuming negligible body forces, azimuthal symmetry and no radial pressure 
gradient reduces Equation 20 to 


du , du 2 , duv , uv 1 dP , ( d 2 u , d 2 u , 1 du \ 

ir + — + — + — = ---dT + \-^T + -^ + -^r) 


( 21 ) 


4.3.2.4 Conservation of Energy 

Full conservation of energy for an incompressible, constant property fluid is 
given by 


_ dT , dT , dT , w dT 

pc {— +u 17 +v — + -^^ 


kV-T+ 2n< 


fr) 2+ 


dy_ 

dr 


+ 


dw 


-=- -^-+v 


d(f> 


+ 


(22) 


dv di u 

dx dr 


+ 


dw 1 du 

dx r d4> J 


+ 


1 dv 
r d<f> 


1 d ( W \ 

+ r TA~) 
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Neglecting viscous dissipation and assuming azimuthal symmetry results in 


5L + „1L + „4Z1 

at ox or 


= a 


a 2 r 

a * 2 


+ 


d 2 r 


Or* 


+ 


1 dT 
r dr 


( 23 ) 


The conservation of mass and momentum equations are solved subject to 
the boundary conditions using the integral formulation developed in 
Appendices B, C, and D. The conservation of energy equation is solved using 
the finite difference numerical formulation developed in Appendix E. The 
computational implementation of the solution approach is given in Section 
4.5.2. 


4.3.3 Boundary Conditions 

The heat pipe boundary conditions can be considered as either external or 
internal. The external conditions refer to the communication of the heat pipe 
external surface with the surrounding environment. The internal conditions 
are of two types: fluid-solid and fluid-fluid. The fluid-solid conditions refer to 
the interaction between the liquid and either the pipe internal surface or the 
wick. Fluid-fluid conditions govern the interaction between the liquid and 
vapor phases. 

4.3.3. 1 External Boundary Conditions 

The heat pipe external boundary can be divided into three convenient 
regions as shown in Figure 15. The evaporator region absorbs heat from the 
environment, the adiabatic region exchanges no heat with the environment, 
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and the condenser rejects heat to the environment. The evaporator surface is 
assumed to be exposed to a known surface heat flux condition, while the 
condenser exchanges heat through radiation 1 to a sink of known 
temperature[34] . The conditions on the pipe external wall are then 


AXIAL RANGE 

RADIAL RANGE 

CONDITION 

0 < x < L e 

r= R w 

, dT , „ 

« a = <1 in(*. 0 

or 

L e < x < L e + L a 

r= K 

o 

II 

b, k. 

fra 

L e + L a < x < L 

r = R w 

,8T _ 

^ dr ^ ra d 


where q in is specified and q rad is to be determined from q rad = eo{ T% — T*) 
assuming environmental sink temperature T e is known. 


1 The condenser could generally exchange heat with the environment 
through conduction, radiation, and convection. However, the study of 
radiation loaded elements is of current research interest for space based 
applications [9]. 
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Long heat pipes with large evaporator and condenser circumferential 
surface areas relative to the surface area of pipe endcaps will normally have 
negligible effects due to conduction in the end caps. With this assumption, the 
boundary conditions are: 


AXIAL RANGE 

RADIAL RANGE 

CONDITION 

o 

II 

R[< r < R w 

ST =0 

dx 

x ~ L 

R, < r < R w 

8T =0 

dx 
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4. 3. 3. 2 Internal Boundary Conditions 

Solid-fluid surfaces correspond to the heat pipe container internal wall 
interface with the working fluid, and the capillary wick interface with the fluid. 
The pipe wall conditions are no slip and matching temperature gradient. The 
conditions are: 


AXIAL RANGE RADIAL RANGE 


CONDITION 


x = 0 R v < r < R, 


u = 0 


v = 0 
dT t 


dx 


= 0 


x = L 


R v < r < R[ 


u= 0 


v = 0 
8T, 


dx 


= 0 


0 < x < L 


r — Ri 


u = 0 


v= 0 

, ST, , ST W 
k ‘ dr k ' dr 


0 < x < L 


r=R„ 


See Section 4. 3. 3. 3 
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4. 3. 3. 3 Liquid-Vapor Interface Conditions 

The liquid-vapor interface of the heat pipe is the most difficult feature to 
analyze. The interface serves as the medium for communication between the 
liquid and vapor along the entire axial length of the pipe. The flow 
characteristics are distinctly different between the respective liquid and vapor 
sides of the interface. The interface is subjected to physical phenomena 
through the action of surface tension that is generally not considered in the 
liquid or vapor flow fields . 2 The general liquid-vapor interface conditions are 
reduced in Appendix A to the following forms used in this analysis. 


KINEMATIC SURFACE CONDITION 

t/ fr = 0 (24) 

where U tv is the velocity of the liquid-vapor interface surface. 


CONSERVATION OF MASS 


Pv V Ov = Pl v 0l 


(25) 


2 The vapor and liquid regions are each assumed to be a single phase such 
that multi-phase phenomena in these bulk flow regions are not 
considered. Entrainment of liquid by shearing action of the vapor, 
condensation in the bulk vapor, and vapor bubble formation in the liquid 
are common examples of departures from the single phase assumption. 
These phenomena are outside the scope of the present analysis. 
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where l' 0v is the vapor radial phase change velocity and P 0/ is the liquid radial 
phase change velocity. In future uses. F 0 will refer to the liquid phase change 
velocity unless for clarity F 0/ is used. 


CONTINUITY OF TANGENTIAL STRESS 

As discussed in Appendix A, the tangential stress condition provides an 
important coupling relation between the liquid and vapor phases. Since the 
vapor phase solution is not included in this study, as a first approximation the 
vapor tangential shear stress will be neglected and replaced with the axial 
velocity no slip condition 

" 0 =° ( 26 ) 


CONTINUITY OF NORMAL STRESS 


K- Pi+ ^v'V'-'ov- 


I'o /) + 2 


dvj 




( 27 ) 


CONTINUITY OF THERMAL FLUX 


dT, 


C »Pv V C>v h fg=- k l--- L 


( 28 ) 


where 
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Cw? _ D p +. D w _ 

The term C w accounts for the effective liquid-vapor interface surface area due 
to the presence of the wick pores of diameter D p and of the screen wires of 
diameter D w . 

SUPPLEMENTAL RELATIONS 

Supplemental relations can be used to derive relations between the liquid 
and vapor at the interface. Using the kinetic theory for interphase mass 
transfer[35] and the Clausius-Clapeyron Equation[36] , the following 
expression is derived in Appendix C for the radial liquid velocity at the 
liquid-vapor interface 

■ =Ai[r v - 3,2 ](7- v -7' J ) 

The interphase heat flux can be expressed by using the continuity equation and 
substituting Equation 29 in Equation 28 

% -■ s/ifc - T s) < 3o > 


-^yJ[nr 3,2 ](7' v - D 


(29) 
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4.4 THE VAPOR PHASE MATHEMATICAL FORMULATION 

Heat pipe dynamics are strongly dependent on the vapor phase so that 
complete solution of the heat pipe problem requires solution of the full system 
of vapor phase governing equations. The intent of this research was to develop 
a computational algorithm for liquid phase dynamics that would be combined 
with a vapor phase algorithm. Development of the vapor phase algorithm 
encountered extreme difficulties such that a vapor code was not available for 
complete checkout of both the liquid and vapor codes[37] . An approximation 
of the vapor phase was developed to permit demonstration of the liquid phase 
code operation. 

The characteristics of vapor flow allow some simplification in developing 
the approximate model. In general the vapor response is orders of magnitude 
faster than the liquid response such that the vapor can possibly be treated as 
quasisteady in a liquid time scale reference [18]. The vapor is also assumed 
to be uniform in a selected control volume. A mass balance on the control 
volume for a single arbitrary time step gives 

Al { n+ 1 ) = -,S t m v (31) 

where M v is the total vapor mass at a given discrete time, 6, is the time step 
size, and m v is the net mass flux due to evaporation and condensation. The 
sign of m v is taken such that evaporation, which has a negative velocity, 
produces mass addition, while condensation produces mass extraction from the 
vapor space. 
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Solution of Equation 31 can be used to fix the vapor state. The solution is 
dependent on the phase change process, which is a function of the heat pipe 
operating limits. This study will consider restriction only due to the sonic limit. 
The solution is then dependent on whether the pipe operation is subsonic or 
sonic limited. 


4.4.1 Subsonic Vapor Operation 

Subsonic operation is assumed to occur such that the vapor is uniform at a 
given time step. The mass balance control volume is then the entire vapor 
space. The net phase change mass transfer can be expressed as 

m v = f p y V^C^L 71 K/bc (32) 

J 0 

Using the continuity equation (Equation 25) and the kinetic energy 
formulation for V 0 (Equation 29) gives 



where T y is the uniform vapor temperature, T s is the liquid-vapor interface 
surface temperature, and R g is the gas constant for the working fluid of 
interest. As a first approximation, all parameters except the temperatures in 
Equation 33 are assumed to be at time level /<">. 

Substituting Equation 33 in Equation 31 and dividing by the fixed volume 
of the vapor space ¥ gives 
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(34) 


p[" + 'W,"V- 


f£;A 

4 -' 




As shown in Appendix F, the Clausius-Clapeyron Equation can be solved for 
the vapor density to give 



All parameters except temperature are again assumed to be from time level 

tin). 

Equations 34 and 35 can be equated to produce a nonlinear expression for 
T v in terms of time level t n) properties and time level liquid-vapor 

interface surface temperatures. The equation can be solved using Newton's 
method. 


4 . 4.2 Sonic Limit Operation 

The sonic limit condition is a restriction on axial mass and heat transfer. 
This restriction implies that the entire vapor space cannot be treated as 
uniform. The vapor space will be approximated as two control volumes. The 
first control volume is the vapor space of the evaporator. The second control 
volume consists of the combined vapor space of the adiabat and condenser 
regions. The two control volumes communicate through the axial vapor 
transport corresponding to sonic limited operation. The control volume mass 
balance for the evaporator is 
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(36} 


M, 


(«+!)_ \fW-b m -6 m 

— M ve °i m ve °t m son 


and for the adiabat-condenser is 


I W<r 1) =A4?-a,m v£ +a> J „„ (37) 

where m son is the sonic limit axial vapor flux. The sonic limit mass flux can 
be derived from the sonic limit heat transfer given by Equation 4. Solving for 
mass transfer of sodium vapor produces 


m 


son 




(38) 


The phase change mass flux for each control volume is found from 
Equation 32 with limits of integration corresponding to the respective control 
volume. The evaporator control volume has phase change mass flux of 


m 


ve 


r L e 

= Pv V 0v C *2 * 


(39) 


or 


m ve 


CwRv J~^ Pvh & 


j 7 ^ — 3/2 

^ e 1 v 1 v 


C L e 

T s dx 


(40) 


Phase change mass flux for the adiabat-condenser control volume is 


m 


VC 


= | Pv V 0v C v?- n 


(41) 


or 
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(42) 


/ 27T 

m \'C — — P\hfg 

Equations 36 and 37 can be solved using Equations 38, 40, and 42 for the 
vapor density during sonic limit conditions for each of the control volumes. 
Using the Clausius-Clapeyron Equation as in the subsonic case produces a 
nonlinear expression for vapor temperature in each control volume that can 
be solved with Newton's method. 


{L 


- L e )r;' A - T; ij2 j L r s dx 


4.5 THE COMPUTATIONAL IMPLEMENTATION 

The governing equations, boundary conditions, property relations, and 
supplemental equations were coded for solution in a Fortran computer 
program. The calculations are performed at discrete grid points in the liquid 
and pipe wall. For clarity, the program First will be presented in a top level 
sense to establish the overall sequence of calculation. The detailed algorithm 
will then be discussed. Equations presented in other sections of the thesis are 
repeated when discussed below for convenience. 

4.5.1 Computational Algorithm Overview 

Overall organization of the program is shown in the flowchart in Figure 17. 
The algorithm begins by setting various parameters and initializing conditions. 
The time stepping section is then entered. Operating limits and property 
dependent parameters are updated. The input heat flux is updated as is the 
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vapor regime. The iteration section is then entered. The first partial time step 
(t*) is solved for liquid and pipe temperatures, vapor temperature, and liquid 
velocities. The same steps are then repeated for the second partial time step 
(/('H-O). Iterations arc repeated on the first and second time split equations 
until temperatures of the liquid surface and container wall have converged to 
a user specified tolerance. After time step convergence the thermophvsical 
properties are updated. Heat transfer limits are determined. Input, output, 
evaporation, and condensation heat transfer rates are calculated. The 
capillary limit heat transfer is calculated and compared to evaporation heat 
transfer. A stop flag is set if the limit is exceeded. A steady state test is 
performed by comparing input and output heat fluxes. Flags are set if steady 
state is reached. The liquid pressure gradient and liquid pressures arc 
evaluated. The capillary radii are determined from vapor pressure, liquid 
pressure, and phase change velocities. Capillary radii are checked for the 
capillary pumping limit. A stop flag is set if the capillary limit is exceeded. 
Time step data are printed, and the program advances to the next time step 
unless a stop flag was issued or time has' expired. If calculations are to be 
terminated, restart data are printed to a restart file which can be used to 
continue the calculations in a later run. 
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4.5.2 Computational Algorithm Details 

The detailed program can be considered in three broad sections: the startup 
section, the iteration section, and the completion section. The startup section 
shown in Figure 18 begins the calculation cycle. Calculations can be 
performed beginning from a time zero condition or in restart mode to continue 
calculations of a previous case. The algorithm begins by setting common and 
dimension areas, and reading input data. Input data consist of the restart . 
option flag, time step size, grid step sizes, vapor space radius, wick pore radius, 
screen wire radius, environment heat sink radiation temperature, and initial 
heat pipe temperature. The restart option flag controls the program 
initialization steps so that for a restart case, the program will initialize with the 
appropriate data from the previous case. Geometry parameters are calculated. 
Time parameters controlling the duration of the calculations are set. Data 
output print control parameters are set. Heat pipe surface radiation 
characteristics are set as are input heat flux control parameters used in 
Equations 44 and 45. Thermophysical properties of the pipe wall material and 
sodium working fluid properties are initialized. Combinations of various 
constants used repeatedly throughout the code are calculated to reduce total 
arithmetic. Peaceman-Rachford ADI parameters used by Equations 134 and 
135 in Appendix E are initialized. Operating limits are initialized, and 
temperature and velocity fields are initialized. Startup data are printed, and 
control passes to the beginning of the time steps in the iteration section. 

The iteration section shown in Figure 19 begins the time step and performs 
iterations between temperatures and velocities until acceptable convergence is 
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achieved. Property dependent operating limits and ADI parameters are 
updated. Input heat flux is updated for the new time step. The input heat flux 
profile and time dependence are user defined functions given by Equations 44 
and 45 

<7 in = y W j 1 - sin k (t + y) } (44) 

t > T 

Q 'm ~ ^max ^ 4 ~^ 

An option is included to smooth the abrupt step change of heat input at the 
evaporator-adiabat boundary. The vapor operating regime is checked using 
Equation 38 



Operation is subsonic if axial heat transfer is below the sonic limit and is sonic 
limited otherwise. 

Iterations are set up by saving data from the just completed time step and 
using these data as the initial guess for the next time step. Iterations are 
performed by first solving for time t* liquid and container temperatures. The 
finite difference Peaceman-Rachford Alternating Direction Implicit method is 
used to solve the energy equation in two partial time steps. The first partial 
time step is solved using Equation 134 from Appendix E 

Afat + A 2 T-j+ A 3 T- j+i = fi, TjV fi 3 Ttfv (134) 
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Implementation details are provided in Appendix F. The method essentially 
involves solution by integrating the energy equation for a partial time step 
(time t*) in the radial direction, and then integrating for the remaining time 
step (time r< n+1 >) in the axial direction. 

The first time split liquid temperatures are then used to find time t* vapor 
temperatures. Depending on whether the vapor is subsonic or sonic limited, 
control passes to the appropriate subroutine to use Newton's method to solve 
the vapor temperature equation. Subsonic vapor temperature is determined 
from Equations 34 and 35 



Sonic limited vapor temperatures are determined from Equation 35 above and 
Equations 36 and 37 

M { Z + ' ] = (36) 

<> -«,< + S, m son (37) 

Based on vapor temperatures and liquid-vapor interface surface temperatures 
the liquid phase change velocity is calculated from Equation 29 
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(29) 


77 /l /s]c r '~ 3,2 ] (r v- T,) 

= *i[77 3,2 ](n- t ,) 

The bulk axial velocity distribution is calculated from the continuity equation 
and the liquid phase change velocity using Equation 84 from Appendix B 

a r 

^') +! Zw + W (S4) 

a L n= 2 

Bulk axial velocity is transformed to a velocity profile using Equation 85 from 
Appendix B 

u* = a + b r* + c r* 2 + d r* 3 (85) 

Finally, the radial velocity distribution is found from Equation 119 of 
Appendix C 



v(x,r) = 


c w y 0 K 


R[ + R v 

R„ \ .*3 


R v r *2 / R \ .3 / R 

VT + ( t + ^T + ( f+ "f 


A? \ *4 *5 
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Calculations for the first time split are completed. 

The second time split is set up by setting the temperatures and velocities for 
the first iteration of the (/<" +1 >) time split equal to the time t* values. The liquid 
and container temperatures, vapor temperatures, and velocities for the second 
partial time step (r< n+1 )) are found as in the first time split with liquid and 
container temperatures found using ADI Equation 135. from Appendix E 
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(135) 


^4 7 1-tJ )+ A 5 7 1j +1)+ /1 6 7 /+tj )_ Vi/-1 + B 5 T iJ + B 6 T iJ+ 1 


The energy equation and the flow dynamics equations are tightly coupled 
at the liquid-vapor interface such that the solutions of each are interdependent. 
In addition, the equations contain nonlinear terms that required 
quasilinearization treatment, which requires iterations. Iterations are 
performed between temperatures and velocities until an acceptable 
convergence is achieved. 

The completion section of the program shown in Figure 20 includes all 
remaining functions to complete the time step calculations and send control to 
the next time step or terminate execution. Thermophysical properties are 
updated based on the temperature conditions using the relationships presented 
in Appendix F. A mass balance is performed to track working fluid inventory. 
The sonic limit heat transfer is determined from Equation 4 
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Capillary limit heat transfer is determined from Equation 8 
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Input, output, evaporation, and condensation heat transfer rates are 
evaluated. A stop flag is issued if the capillary limit is exceeded. Input heat 
transfer is compared with output heat transfer, and a steady state flag is issued 
if they are approximately equal. 



Figure 20: Algorithm Completion Flowchart 

Liquid pressures are found from Equation 126 in Appendix D from the 
axial momentum integral equation 
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The integration is performed by assuming liquid pressure and vapor pressure 
are equal at the end of the condenser ( x = L). Details of the numerical 
implementation are provided in Appendix D. The liquid-vapor interface radial 
momentum condition given by Equation 27 can be solved for the capillary 
radii of curvature to give 



The total liquid pressure drop is compared to the maximum possible capillary 
pumping capability of the wick given by Equation 10 

(10) 

A stop flag is activated if the limit is exceeded. Data from the time step are 
printed, and the next time step is initiated unless a stop flag has been activated 
or time for the case has expired. Data from the final time step are dumped to 
a restart file, and execution is terminated. 
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Chapter V 

DISCUSSION OF RESULTS 


“It is of highest importance in the art of detection to be able to recognise 
out of a number of facts which are incidental and which vital. Otherwise 
your energy’ and attention must be dissipated instead of being 
concentrated." 


Sherlock Holmes in The Adventure of the Reigate Squire 


Capabilities of the computational algorithm were tested by running two 
heat pipe startup cases. The heat pipe dimensions were chosen to be consistent 
with the experimental device described by Merrigan et al. up to the limit of the 
uniform grid used in this analysis [19][20]. The experimental heat pipe in 
that study used lithium as the working fluid and an annular wick 
configuration. The annular region for liquid flow was formed between the pipe 
interior wall and a porous concentric tube constructed of 7.25 layers of pressed 
screen. The high temperatures involved required the container to be 
constructed from the refractory alloy molybdenum. The heat pipe geometry 
parameters used in this analytical study compared to the actual parameters are 
(see Figures 15 and 16): 
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Parameter 

Actual 

Approximate 

Heat Pipe Internal Length (L) 

4.0 

4.0 m 

Evaporator Length ( L e ) 

0.4 m 

0.4 m 

Condenser Length ( L c ) 

3.0 m 

3.0 m 

External Diameter ( D w ) 

1.90 cm 

1.886 cm 

Internal Diameter (Dj) 

1.60 cm 

1.598 cm 

Vapor Space Diameter ( D v ) 

1.49 cm 

1.490 cm 

Effective Wick Pore Diameter ( D p ) 

i 53 nm 

53 ^m 


where the effective wick pore diameter was experimentally determined from 
surface tension measurements. 

An axial grid step size of h x = 0.2 m and a radial grid step size of 
h r = i.i x 10~ 4 m were used. The computational grid consisted of 21 axial grid 
steps and 20 radial grid steps. The grid mapping to the heat pipe configuration 
is: 


Radial Liquid Space \ <j<6 

Radial Pipe Wall 6<y<20 

Axial Evaporator Region 1 < / < 3 

Axial Adiabatic Region 3 < i < 6 

Axial Condenser Region 6 < i < 21 


Numerical experimentation showed that the code is highly sensitive to time 
step size. Too large a time step produced an inconsistent temperature at the 


84 


liquid-vapor interface. The temperature was inconsistent in that the 
temperature gradient at the surface changed sign relative to the gradient of the 
interior liquid. For example, in the evaporator energy is conducted from the 
pipe wall to the liquid-vapor interface. The liquid-vapor surface temperature 
should be the lowest temperature in the radial conduction path in order for the 
energy to reach the surface. Too large a time step will produce a physically 
questionable surface temperature that is greater than the temperature of the 
interior liquid. Too large a time step was also found to cause oscillations in the 
pressure calculations. 

The time step size was selected by running a series of calculations for three 
hundred time steps each with gradually decreasing step sizes. A time step size 

<5 f =2x 10~ 3 second was found to eliminate inconsistent liquid-vapor 
interface temperatures, but did not completely remove pressure oscillations. 
A time step size of 5 , = 1 x 10 -3 second removed the pressure oscillations for 
the three hundred time steps calculation. This time step was used in the 
calculations. 


5.1 SLOW STARTUP TEST 

A startup test was devised to supply input energy to the evaporator such 
that the heat pipe would warm up without reaching the sonic limit. This test 
was to check overall program implementation. Input heat flux was supplied 
according to the function 
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t < T 


^ in 



1 — sin 




(44) 


t > T 


Qin ~~ ?max 


(45) 


where q miX is the maximum input heat flux and t is the time stretching 
parameter that controls how quickly q m reaches q max - In this test, 
n = 1 x 10 6 IVI m 2 and t = 300 seconds corresponding to a five minute 

H max * 

period for the input heat flux to reach the maximum. 

Calculations were performed for 114,193 time steps, or a period slightly 
greater than 1 14 seconds. The run was terminated at this time since the liquid 
pressure drop exceeded the capillary pumping limit. The liquid pressure drop 
and the capillary pumping limit as functions of time are shown in Figure 21. 
Analysis of the results shows that the capillary pumping limit was exceeded 
due to a large jump in pressure drop at a single time step. Figure 21 also 
shows that the pressure calculations had been experiencing oscillations with 
usually small amplitudes of approximately 2.5 Pa. The oscillations began after 
the 300 time steps used to verify the adequacy of the time step size of 
1 x 10 ~ 3 second. The liquid pressure drop curve is drawn through the centers 
of the oscillation ranges in Figure 21. The oscillations become very large at 
t _ 90 seconds, but begin to decay until the large jump occurs at approximately 

1 14 seconds. 
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The pressure oscillations are presented in detail in Figure 22. The liquid 
pressure drop is plotted for eleven consecutive time steps for 
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19,990 < k < 20,000, corresponding to the time around time /= 20 seconds in 
Figure 21. The oscillations are of essentially constant amplitude and well 
behaved until time t > 80 seconds as shown in Figure 21. 
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Figure 22: Slow Startup Liquid Pressure Drop Oscillations 


Power transfer as a function of time is shown in Figure 23. While the input 
power briefly exceeds the sonic limit line, the phase change heat transfer 
remains below the sonic limit. At approximately 103 seconds into the 
transient, the sonic limit exceeds the capillary limit so that the capillary limit 
becomes the controlling heat transport limit. At approximately 105 seconds, 
the unlikely result is produced that the phase change heat transfer exceeds the 
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input heat transfer. With the exception of this anomaly, the heat transfer 
curves are well behaved, smooth functions of time. In addition, the phase 
change heat transfer is well below the capillary limit heat transfer at the point 
when liquid pressure drop jumps to exceed the capillary pumping limit. The 
cause of the liquid pressure drop jump is not apparent from the heat transfer 
curves. 

The vapor temperature as a function of time is shown in Figure 24. The 
vapor temperature is also a smooth, although rapidly increasing, function of 
time. The cause of the jump in liquid pressure drop is also not apparent from 
vapor temperature. 

The pressure drop calculation can be seen in Figure 21 to produce 
oscillatory behavior with large sporadic oscillations for time greater than 80 
seconds. The pressure calculations are a sensitive function of time step size. 
Numerical experimentation also showed that the pressure calculations are 
sensitive to property updates and changes in the liquid-vapor interface 
boundary conditions corresponding to the change between subsonic and sonic 
limited operation. Sensitivity to property updates was identified by varying 
the frequency at which updated properties are calculated. 

There are many potential causes for the pressure oscillations. The pressure 
oscillations may be physical, but this possibility is remote considering the 
numerical sensitivities of the method. The time step size of 1 x 10~ 3 second 
may have been too large. If this is the case, a smaller time step may not be 
acceptable considering the cost of computations for a transient process that 
occurs over a period of tens of minutes. The oscillations may have resulted 
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from the derivation in Appendix B of the axial velocity profile in which the 
shape is determined by neglecting the temporal terms in the momentum 
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equation. The pressure calculations are also based on the convenient 
assumption that the minimum liquid pressure occurs at the end of the 
condenser, although this assumption is not generally valid for the wide range 
of operating conditions encountered during a startup transient. 
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5.2 SONIC LIMIT STARTUP TEST 

A startup test was performed to test algorithm calculation of sonic limited 
operation. Using again r = 300 seconds, the maximum input heat flux was 
increased to q m ^ = 1.5 x 10 6 Wjm 2 to assure the phase change heat transfer 
would reach the sonic limit. Heat transfer results are shown in Figure 25. 

Input power exceeds the sonic limit at 20 seconds, and remains above the 
sonic limit through the duration of the test of 130 seconds. The evaporation 
and condensation heat fluxes are essentially equal as the sonic limit is reached 
at approximately 40 seconds. Evaporation and condensation follow the sonic 
limit curve until 80 seconds. After 80 seconds, condensation heat flux 
increases above the sonic limit, while evaporation heat flux continues to follow 
the sonic limit curve. This unlikely result is a result of the inadequacies of the 
simple vapor model used to perform calculations. The vapor model 
unintentionally forces a high condensation rate, which forces the condenser 
temperature to increase rapidly due to the high thermal resistance of the 
radiation boundary condition on the condenser. This temperature effect is 
shown in Figure 26. The evaporator vapor temperature smoothly and 
gradually increases with time, while the condenser vapor temperature rapidly 
increases. The condenser temperature eventually increases sufficiently so that 
radiation heat transfer from the condenser becomes very efficient. At 
approximately 93 seconds, the output heat transfer exceeds the evaporation 
heat transfer, which is another unlikely result. Calculations are terminated at 
1 30 seconds due to the condensation heat transfer exceeding the capillary heat 
transfer limit. While the capillary limit is not strictly a limit on condensation, 
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Vapor Temperature, 
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Chapter VI 

CONCLUSIONS AND RECOMMENDATIONS 


The heat pipe is a very complex device such that performing analysis of heat 
pipe transients is a difficult venture. The difficulty is due to the coupled liquid 
and vapor dynamics, the phase change process, the coupled heat transfer 
problem with nonlinear boundary conditions, the physical geometry, and the 
transient time requirement. This study has presented the full system of 
governing equations, including boundary conditions, required to solve the 
liquid phase heat pipe problem. A simplified solution was formulated by using 
certain assumptions and integrating the liquid dynamics equations. 

The resulting formulation identified a key requirement for future research 
in solving the heat pipe problem. This study used a kinetic theory approach 
to model the phase change process as expressed by Equation 29. The model 
consists -of a large coefficient multiplying a very small temperature difference. 
The model requires a very small time step for stability even with the 
simplifications used in this analysis. The large coefficient also transforms 
temperature differences that are otherwise beyond machine accuracy into 
computationally significant terms. Future research should investigate an 
alternative model for the phase change process. 

Difficulty with geometry was encountered due to the widely disparate 
length scales in the different coordinate directions. The radial dimension 
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across the liquid is much smaller than the radial dimension across the pipe 
wall, which in turn is much smaller than dimensions in the axial direction. 
One approach to address this problem is to assume radial gradients are 
negligible. The radial pressure gradient was neglected in this study. This 
study also found small radial temperature gradients indicating that the radial 
thermal resistance is small. The drawback of this approach is the loss of 
fidelity with the true physics of the problem. A preferable approach may be 
to use a variable mesh grid with a to-be-determined simplified treatment of 
radial gradients. The variable mesh grid allows fine resolution in important 
regions with coarse resolution in the other regions. The simplified treatment 
of radial gradients would take advantage of the small radial gradients. 

The long duration of transients presents a competing concern with model 
complexity. An overly complex model with high fidelity may be too 
computationally expensive for practical calculations of transients. 

This study has also shown that the transient heat pipe problem is not 
amenable to a straightforward simplification such as is used for the vapor 
phase. The solution of the heat pipe transient problem requires a full solution 
of the liquid and vapor phases. This elusive solution is left to future efforts. 
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Appendix A 

LIQUID-VAPOR INTERFACE CONDITIONS 

The liquid-vapor interface is a complex boundary that allows 
communication between the liquid flow and vapor flow along the entire length 
of the heat pipe. The complexity of the interface physics can be conceptualized 
as shown in Figure 27. Solution for the flow Field on one side of the boundary 
is dependent on the flow Field on the other side of the boundary. The location 
and conFiguration of the interface surface fluctuates with operating conditions. 
The phase change process, including direction of phase change, is a function 
of conditions on both sides of the boundary. The interface can be seen to 
involve exchanges of mass, momentum, and energy. 
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In general two transition relations are required to define the behavior of 
transportable quantities at material surfaces[38] . One relation is required to 
represent continuity of intensity of the appropriate quantity across the surface 
while the other relation represents continuity of the normal component of the 
flux vector across the surface. 

The principles of continuity of intensity and continuity of normal flux are 
not generally exact for non-equilibrium processes. Net evaporation and 
condensation, which occur simultaneously in different regions of a heat pipe, 
are by definition non-equilibrium process. There is actually no rigorous 
justification for requiring, for example, matching liquid and vapor 
temperatures at the interface. Similarly, the requirement for matching 
velocities (the so-called no-slip condition) at fluid-fluid boundaries such as the 
liquid- vapor interface in a heat pipe is not as certain as is the no-slip condition 
at a solid-fluid interface. Errors introduced by non-equilibrium conditions are 
dependent on the severity of the departure from equilibrium. In many cases, 
as is assumed for the heat pipe, departures from local equilibrium are small 
such that non-equilibrium errors can be neglected. The equilibrium condition 
can further be justified based on the action of diffusive transport processes. 
Diffusive forces resist non-equilibrium conditions in proportion to the severity 
of the departure from equilibrium. Any non-equilibrium condition will tend 
to be forced back to equilibrium such that non-equilibrium errors will be small. 

The full set of conditions for free liquid-vapor interfaces is given by 
Prosperetti[39] . The interfacial region separating liquid and vapor is assumed 
to be a massless two dimensional surface across which transport quantities 
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undergo an abrupt change. The actual physical situation involves a transition 
layer between liquid and vapor across which there is a rapid but continuous 
change of transport quantities. The jump conditions approximate the change 
in flow quantities in the limit of an infinitely thin transition layer. 
Inconsistencies with the jump conditions can be reconciled by assuming certain 
surface sources for momentum and energy related to surface tension. The full 
set of interface conditions are reduced below using the assumptions 
appropriate for this analysis. 


A.l KINEMATIC SURFACE CONDITION 

The problem involving a surface free to evolve depending on local 
conditions requires determination of the shape and position of the surface with 
respect to time. The kinematic boundary condition is used for this purpose. 
Let the liquid-vapor interface surface be governed with respect to time by the 
function 5 such that the shape of the surface is given by 

5 (*, f)=0 

The kinematic boundary condition is then 

— + U h ;VS = 0 on 5=0 (46) 

dt n 

where U lv is the velocity of the liquid vapor interface surface. Systems far 
removed from equilibrium may have violent evaporation or condensation such 
that the surface velocity is of the same order of magnitude as the flow field 
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velocities. Conversely, surfaces in systems close to equilibrium such as is 
assumed for the heat pipe will have negligible velocity. When the interface is 
stationary or has negligible velocity, the condition reduces to the usual 
requirement that the normal component of velocity vanish 

U h :n = 0 (47) 

A. 2 CONSERV ATION OF MASS 

The general conservation condition is 

Pl(U 0! - U^'n = p v (U Qv - vT (48) 

where U 0l and U 0v are the surface velocities of the liquid and vapor 
respectively, p l and p v are the liquid and vapor densities, and n is the unit 
normal to the liquid-vapor interface surface. 

Neglecting the velocity of the liquid-vapor interface gives 

p / U or n=p v U Qv ‘n (49) 

Further refinement is achieved by assuming that there is no axial slip at the 
liquid- vapor interface relative to the wick as discussed in Section A.3, and 
there is no azimuthal flow. Additional simplification is possible by assuming 
that t x n r Pt r n x such that t x n r cz 1. This assumption is equivalent to 
approximating the liquid-vapor interface as a cylindrical surface such that 
velocities at the surface are in the radial direction only. Heat pipe regions with 
a flooded wick such as portions of the condenser where there is little difference 
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between liquid and vapor pressure would satisfy this approximation closely, 
while the regions where liquid recedes into the wick would only approximately 
satisfy this requirement. Conservation of mass reduces to 

P/l'oZ-Pv*' Ov ( 5 °) 

where V ol and V 0v are the liquid and vapor radial velocities at the surface. 


A.3 CONTINUITY OF TANGENTIAL STRESS 

The general condition is 

n x(Z/ 0v - P fr )[p v (i7ov- Uk)-n~\-n x[(T v -T/)*n] = « x (51) 

where t v and t, are the viscous stresses and V r is the surface operator. 

The right hand side gives the tangential stress contribution of surface 
tension gradients. The surface tension of a pure substance is generally a 
function of temperature. Convective force due to surface tension gradients 
resulting from local temperature gradients is known as Marangoni Convection. 
Local temperature gradients in the tiny screen pores at the liquid-vapor 
interface are assumed small such that Marangoni Convection effects can be 
neglected so that the right hand is zero. 

The first term on the left hand side is also identically zero due to the no-slip 
condition at the wick discussed later, nx(U v - £/,)= U 0v — U ol = 0 . The 
remaining terms of Equation 51 give an expression for tangential stress at the 

interface due to viscous forces 
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(52) 


n x [(t v - T/) • n ] = 0 

The total stress vector P has tensor components 

P i =Pn i -o ik n k (53) 

The tangential stress is formed by the scalar product of P and the unit tangent 
vector to the liquid-vapor interface surface t 

t -P = t i P i = p t . n . - t . o ik n k (54) 

where tensor summation notation is used. 

Note that t ( n t ■, = t *n — 0 since the unit normal and tangential vectors are 
orthogonal. Also for i = k Equation 54 becomes t k P k = - o kk t k n k , but 

n k ~ t " n = 0 so that the i = k terms vanish in writing out the components 
of tangential stress 


' • P ~ ~ *x a xr n r ~ ?x °x4> n 4> rx n x ~^ r<f> n cj> ^ x n x a<f>r n r (55] 

The stress elements are given by Landau and Lifschitz as[40] 


_ _ _ ( du dv \ 

xr rx ~ "( dr + ~d7) 

n __ _ . / dw 1 du\ 

X <P <t>x py dx + r 


„ _ _ _ 1 dv , dw w 
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Invoking the assumption that there are no azimuthal flows or gradients for the 
assumed two dimensional flow gives 


o 


xr 




a x4> 


= 0 


0 *r<f> 


= 0 


Tangential viscous stress is now 


t ‘P — a xr (t x n r + tf n x) 


(56) 


Bv using the cylindrical surface approximation, tangential viscous stress is 


given by 


/ -P 



+ 4 ^) 

ox J 


(57) 


The continuity of tangential stress condition finally reduces to 



Equation 58 provides an important relation for coupling liquid and vapor 
dynamics, particularly for determining liquid flow reversal and entrainment 
due to vapor shear forces. However, as discussed in Section 4.4, a vapor model 
was not available for coupling to the liquid model. In order to perform 
computations, viscous stress in the vapor is neglected. This assumption is 
strictly valid only for heat pipe operation at relatively high vapor pressure with 
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vapor speeds that are low such that entrainment concerns can be neglected. 
Invoking the assumption of no viscous stress in the vapor reduces the 
tangential stress condition to the simple result 


Pv 


SJL + Il.) =0 

Or dx A- 


(59) 


Equation 59 does not in general provide a useful boundary condition for 
solving the liquid phase dynamics. The tangential stress condition is replaced 
by the no slip condition at the liquid-vapor interface. This approximation is 
equivalent to assuming that due to the fineness of pore dimensions involved 
(e.g., 30 - 50 /rm pore diameters [9] ), the no slip condition on screen wires is 
assumed to dominate over the free slip condition between screen wires. The 
boundary condition is 


“0 = 0 


(60) 


A.4 CONTINUITY OF NORMAL STRESS 

The general form of the normal stress condition is 

Pv [(^0v — Uh) ' n ][(^0v ~ £/ 0 /) ’ n J + Eh */*]„—[« * p]/ = oV • n (61 ) 

where Vv? is the surface divergence, which is related to meniscus radii of 
curvature and R 2 by 

V-"=-F-+-F- (62) 
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The first term on the left hand side in component form neglecting the surface 
velocity is 


Pv[^0v ' n ][( L '0v Voi) ' n ] _ (63) 

pJyU x n x + \\n r + w\/i^)[{u v - uj)n x + (v v - v,)n r + (w v - »v/>ty] 

Applying the no-slip condition at the wick and no azimuthal velocity gi\'es 

Pvft/ov * n ][(?0v -%)*«] = PvK^x + V*rXK v v - V /K (64 ) 

Assuming that the liquid-vapor interface is approximately cylindrical gives 

p v [ Ufo'it ] [( U 0v - UqiYh ] = P v v( V 0v ~ Vu) ^ 6 5 ^ 

The second term on the left hand side of Equation 61 gives the normal stress 
contribution from the pressure and viscous forces 


n -P = n t Pj = rij P n , - n t a ik n k = P- rij a ik n k 


(66) 


As with tangential stress, o = 0 for / — or k — (f> and o xr — o rx so the normal 
stress components are 


2 2 l 

n \ °ik n k = °xx n x + G rr n r + a U n 4> + 2 ° *r n x n r 


(67) 


where 


a xx 2 P 


du 

dx 


( 68 ) 


~ Ov 

a rr = 2 P-^~ 
Or 


(69) 
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~ / 1 dw . v \ ~ v 

^ = + T j = ^ T 


(70) 


Using the cylindrical interface approximation ( n x = 0) gives 


n-P = P-2fi— 

dr 


(71) 


The continuity of normal stress is finally 


n ~ + PyJ 7 0l ( I'ov - v 0 j) + 2 


dv/ dv v 

Pi— Py— 

Or cr 


,1^1 

=< 7 [ t ; + t - 


( 72 ) 


A.5 ENERGY BOUNDARY CONDITION 

The continuity of thermal energy intensity requires matching temperatures 
at the liquid-vapor interface 
T, 

The continuity of thermal flux condition is 

Pv[(Vov - - */+ y[(^0v- ^) 2 - (^0/- ^) 2 ]} 

+ [(^o/“ U h>y r i]' n - [(^0v - ^)* T v]' w + (?v - ?/)*« = (73) 

1 [ J 57 'dr + + x n )] 

The first term is the energy required for phase change plus the work to adjust 
densities and pressures during phase change plus the energy to accelerate from 
liquid velocity to vapor velocity. The second term is viscous dissipation, taken 
to be negligible at the low velocities in this analysis. The third term is the 
thermal energy conduction. The right hand side is the temperature times the 
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increase of entropy per unit area. The right hand side can be simplified as 
follows. Surface tension is first considered a function of temperature only for 
a pure substance 3 so that 


D do _ d 2 o DT_ ( 74 \ 

Dt dT dT^ 

But cT-ojdT 1 is zero, and if the surface has negligible velocity, the terms in 
parenthesis on the right hand side can be taken as zero so that the entire right 
hand side vanishes. The energy boundary condition reduces to 

p v ( U 0v -n )|/i v - h r + ±[ul - f7o/]} + (?v -?/)■"= 0 ( 75 > 

but 

• n = u x n x + v* r + nyty - 

hy — h[— hjg 

Uq v ~ Ul = Uy + Vfo + Wy - uf - Vqi - wf = f’ov ~ v ll 

{qy-qi)'n=q v .ni-qi.n i 


3 The surface tension as a function of temperature is given in Appendix F 
as 


o = 0.23402- 1 x 10 4 T 
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Assuming heat conduction by vapor is negligible relative to the latent heat of 
vaporization gives 


{q v - q/)-n = -q r n i 


cT, 


= k, 


Ox 

cT, 


Or 


l cTj 57) 

- n I+ —„ r+ —^ 


(76) 


where the assumption is made that radial heat conduction in the liquid is 
dominant over heat conduction in other directions. The energy boundary 
condition reduces to 


^wP v ^ Ov 


1 T ”1 57" ) 

h fg + y < 1 o, - v oi>\ = - k i-rjr 


(77) 


where C w has been added to account for the actual liquid-vapor surface area. 
The kinetic energy contribution at low velocities is negligible relative to latent 
heat of vaporization so that the energy boundary condition becomes finally 


57 ) 

C-wPv I ’th-hfg = ~ — 


(78) 


or using the continuity equation to substitute for the vapor parameters, 


57 ) 

C *Pl l 0 h fg = - k[— 


(79) 
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Appendix B 

AXIAL VELOCITY DETERMINATION 
B.l BULK AXIAL VELOCITY 

The bulk axial velocity V can be found as a function of axial position from 
a control volume mass balance. Consider as the control volume the cell shown 
in Figure 28. Assume the radial interface velocity is uniform along the region 
from i - Vi to i + V» and is equal to known V 0 {i) . Assuming constant liquid 
density in the cell, the bulk axial velocity across cell boundary / + 1 is 

U(i + 1 ) = £ 7(0 + 4 - 4 ^[ V 0 (i) + V 0 (i + 1 )] ( 80 ) 

~ ™ a 

where the cell liquid-vapor interface area A h and the annulus cross-section 
area A a are given by 

Afy = 2nR v h x C w 

A„ = Rh 

The term C w accounts for wick porosity. A solid boundary exists at i = 1 so 
that £7(1)= 0 and the bulk axial flow across grid line i = 2 is 

y(2)=-t--^-[l / oO)+ v^2)\ (81) 

Z A a 

The bulk flow across / = 3 is 
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i/(3) - U{ 2) + -y —y[ v^2) + F 0 (3)] (82) 

^ H a 

Substitution of Equation 81 into 82 gives 

^ (3 )=y-^[^l ) +2r (l (2)+ Fo(3)] (83) 

The bulk flow 7 across interior grid i can then be found from the recursion 
formula 



i — i 

W)+2 '£v<P)+ V^i) 


(S4) 


A solid boundary exists at / = /V so that U(N) = 0. 
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B.2 AXIAL VELOCITY DISTRIBUTION 

The objective is to define the axial velocity in terms of the bulk axial velocity 
U{x) and the radial position r. The velocity profile is assumed to be well 
approximated by a cubic relation 

u* = a + b r* + c r* 2 4- d r* 3 (85) 

where 

* — E 

u ~ U 

* r ~ R v 

r =~T- 

5 = Rj — R v 

Apply the no slip condition at the liquid-vapor interface 

at r = R v , r* = 0, 

u = u =0 
so o = 0 

The three unknowns b , c, and d remain. The no slip condition at the 
liquid-pipe wall interface is applied 

at r = Rj, r* = 1, 
u = u* = 0 

so b + c + d = 0 (86) 
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Equation 86 provides one of the equations required to find the three 
unknowns. Two additional relations are required. One of the required 
relations can be derived by taking advantage of the assumption that there is 
no radial pressure gradient. The axial momentum equation can be solved at 
the liquid-vapor interface for the pressure gradient. This expression for 
pressure gradient will explicitly contain the influence of the phase change 
velocity. This pressure gradient can then be equated to the pressure gradient 
resulting from integrating the axial momentum equation in the radial direction. 
The integral formulation will also contain explicit terms for the liquid-vapor 
interface phase change condition. The resulting expression is the axial 
momentum compatibility condition. This approach is intended to preserve the 
influence of the phase change process on the velocity profile. The time 
dependence of the velocity profile is assumed to be negligible as a first 
approximation. 

The steady, two dimensional momentum equation is 


u 


du 

dx 



1 dP 
P dx 


+ v 



+ 



JL — 

r dr 


Solving at the liquid-vapor interface for pressure gradient gives 


1 dP 
P dx 




_} ^j0_ du 

R v v dr 


K 


( 87 ) 


( 88 ) 


Now the axial momentum equation is expressed in integral form 
corresponding to the heat pipe system of this study. The integration will 
include the temporal terms for general validity since the integral formulation 
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is also used for the pressure determination in Appendix D. The temporal 
terms will then be dropped for the current purpose of determining the axial 
velocity profile. The conservative form of the axial momentum equation for 
two dimensional flow is 

du du " duv _uv_ 1 dP 

dt dx dr r P dx 


~2 ~,2 

d u 6 u + 


du 


dx‘ 


dr- 


dr 


(89) 


Integrate between the heat pipe limits R v and R, 


-Lf* 


urdr + 

cx J 


*i 


u z rdr + 




‘ R i ~ rR, 

rdr + f uvdr = 

dr J R, 


J R V 


1 dP 

P dx 


cRt 

rdr + v 

a 2 

r R i 

urdr + 

r*' 2 

— —rdr + 

r R i * 

du dr 

J *v 

dx 2 ' 

U J 

dr 2 

R v 

*, dr 


(90) 


Integration by parts and evaluation of integrals produces 



Assuming now that the temporal terms can be neglected in determining the 
velocity profile, Equation 91 can be rearranged to give 


1 dP 
P dx 
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Equating expressions 88 and 92 gives after rearranging 



Equation 93 is now manipulated by transforming to u* and r*, substituting 
u* = b r‘ + c r* 1 2 + d r* 3 , and rearranging. The result is 


+ *y 2 

Wi+ K) 


R/+ R , 


R;(b + 2c+3d) 


2 3 




3 S 


H* 2 i+i- : fl + M4+44M + Ml-!-+4 


2 5 


J- I 

5 4 s /J v ax 2 

Vi . 2 M 


5 5 


V 6 5 6/ \ 7 33/ \ 8 Id! v <?jc 


After an enormous amount of algebra and using Equation 86, the following is 
obtained 

T l4 b+ T l5 b 2 + T l6 d + T l7 d 2 + T ]S bd= 0 (95) 

where 



1 1 *v~| <S 3 dhj_ 

12 6 ^ J U dx 2 



1 1 1 2<$ 3 du_ 

60 30 3 J v dx 
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T l6=- 2R l~ R v + 



12 


*v~ 

d 2 d 2 U 

s _ 

U dx 2 


T \7 ~ 


1 1 ~1 2<5 3 dU 

168 + 105 6 v dx 


^18== 


2 + 1 *v 


105 


30 5 


2<5 3 dU 
v dx 


There are now two equations, 86 and 95, for the three unknowns b , 
The third required relation is from the definition of bulk axial flow, 


V = 


rR t 

ulnrdr 



C R 1 

2nrdr 


— 2 f*/ 

U = - - - urdr 

Rf - r; j r v 

Transforming to u * and r* in Equation 97 gives 

*' 2 ~ *1 = S f 'u'r’dr' + R v f \’dr' 

2d J 0 Jo 

Substitution for u* gives 



Rearrangement gives 


■, and d. 


(96) 


(97) 


(98) 


( 99 ) 
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b = 6 + — 
5 


3 - 


R, 


R l + R v 


(100) 


Substitution of Equation 100 in Equation 95 and lengthy algebra result in a 
quadratic equation for d 


S j d + S 2 d "b S 3 = 0 


where 


( 101 ) 




25 


R r 


R I + R y 


1 , 1 R v 


+ 


60 30 <5 


+ 


+ 


1 


H, 


168 105 3 


1 

r 3 - ^ ] 

2 1 R v ~ 

X 

“ « 2 c„^„ 


5 

R/+ R v _ 

_ 105 30 5 

J 

R/+ Ry 

V 


_1_ 

3 

~ r ± 

! -*r 

n 

i n \ 

5 j 

R i + R v _ 


2 Rv 


| 

+ < 

f 1 1 R v 

i 

e 

R v ] 

1 1 

1 

I 20 12 b 

s. 

5 

R l + J 

12 6 b 


R, 


-"-Rf-Rp 


I 


Rj+ R v 

K. 


+ 


12 


1 R v 


b 3 d 2 u 

U ex' 1 


105 5 b 


5 \ R/+R v 


R„ 


5 + 5 <5 


Mcjiv r 0 


R[ + R v 




+ 


dx* 

6+12 4 ^ 2 ^ - 3(^ 2 - tf 2 ) 


<5/5 R[ + R v 
Solution of Equation 101 for d is given by 
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( 102 ) 


d= - S 2 ± V 5 2 _ 4S l s i ] 

The only remaining issue is selection of the proper + operation in Equation 
102. By assuming backflow generally does not occur, the requirement can be 
imposed that the sign of u(x,r) and (J{x) must be consistent at every axial 
location x. Analysis of Equation 102 is not straightforward due to the 
competing signs of the many terms. Numerical experiments were performed 
to study the effect of the operation on the profile. The experiments considered 
using a single operation for the entire pipe, as well as operations dependent on 
the sign of the phase change velocity. These tests showed that for reasonable 
velocities, the only operation that produced no backflow was subtraction. 
Using subtraction in Equation 102 allows the coefficients b , c, and d to be 
evaluated from the known geometry parameters and flow conditions U(x) and 
F 0 (x) using 

*= 6 + t [ 3 - ] rf 

c = —b — d 


(103) 

(104) 

(105) 
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Appendix C 

RADIAL VELOCITY DETERMINATION 
C.l RADIAL PHASE CHANGE VELOCITY 


The objective is to determine the liquid phase change velocity V 0 . A 
straightforward model of phase change mass flux is given by Schrage as the 
simple theory of interphase mass transfer of a pure substance [35] 


Pv ] 0v = 



( 106 ) 


where P s and T s are values associated with the liquid-vapor interface surface. 
The sign of velocity as applied to the coordinate system used in this analysis 
has the convention 


L 0 > 0 Condensation Occurs 

Cq < 0 Evaporation Occurs 

Conservation of mass at the interface from Appendix A is given as 

Pv V i Ov - Pi v 0l 

Assume the vapor is saturated at the vapor temperature 

?v = ^sat I T 

* V 

Assume the temperature jump 5T lv from the surface to the vapor is small 
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T,= T t +6T k 
ST k . « T v 


so that 



1 



Equation 106 can now be written as 


Pf 0 = 




(107) 


Since the temperature jump 5T h . is small, the surface pressure can be 
approximated as 


p ~ p 
r s — r sat 


T 

* V 


dP 


sat 


dT 


(T v - T s ) 
T 

1 V 


so that Equation 107 becomes 



1 dP sat 

JrX dT 


(T v - 

T 

1 V 



(108) 


The slope of a phase curve in a pressure - temperature diagram can be 
related to the latent heat and to the volume discontinuity by the 
Clausius-Clapeyron Equation [36]. In terms of molar entropy 5 and molar 
volume V, 


dP AS 
dT AV 


009) 
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Note that AS and A V are discontinuities in molar entropy and molar volume 
associated with phase change. The latent heat of vaporization is defined as 

^T&S 


So that Equation 109 can be written as the Clausius-Clapeyron Equation 


dP h fg 
dT TAV 


( 110 ) 


or 


dP 

dT 



( 111 ) 


where the right hand side is evaluated at an appropriate temperature, taken in 
this case to be the originating temperature. Substituting the 
Clausius-Clapeyron Equation in Equation 108 produces 


/ M 

1 

1 

/ 2 nR u 

" J 


T v 3/2 (7 v - T s ) 


0r 

*'o-*|77 3,2 d'v- T s ) 

where the assumption has been used that p v <$ p, . 


( 112 ) 


(113) 
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C.2 RADIAL VELOCITY DISTRIBUTION 

The objective is to determine the radial velocity as a function of position. 
The phase change velocity V 0 {x) is assumed to be known from the kinetic 
formulation of Equation 113. The radial velocity distribution can then be 
determined by integrating the continuity equation in the radial direction 
between the limits R v and r. 


rr 


rr 


R v 


cu 

dx 


2nd dd + 




cv 

dr' 


2nd dr' + 


J> 

J R V 


2nd dr' = 0 


(114) 


or 


v(r) = 


£wVogv 

r 



(115) 


The screen permeability factor C w accounts for the liquid-vapor surface area 
reduction due to the presence of the screen. Evaluation of Equation 115 at 
r = R v produces the discontinuous result v(R v ) = Q v o * s0 thal Equation 115 
is not strictly valid at r = R v . 

Evaluation at r — R t , where v (/?/)= 0 due to the impermeable wall 
boundary condition, produces 


V 


o- 


__J d_ 

C W R V dx 


C R l 

uddd 

J R V 


(116) 


A useful relation can be derived by using the expression for bulk axial velocity 
developed in Appendix B 
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r Rl 

- — - ur'dr' 

_ c > 2 Jj? 


Rf - r; j r 


Using U in Equation 1 16 gives 


v 0 = 


1 O 

Rf - r; 


V ~ A V dU 

2C W R V dx 


Returning to Equation 1 15 and transforming the integral to u * and r* produces 


v(r) = 


C w V 0 R v S dU 


— -|<5 J u*r*'dr*’ + /? V J u*dr*' j> 


(118) 


Substituting u * = br‘ + cr * 2 4- dr * 3 allows evaluation of the integrals 


rr* *3 *4 »5 

u*r*'dr*’ = b^—+cl—+dl— 
J n 3 4 5 


rr* *2 *3 *4 

u *dr*' = b ■!— + cl— + rf— 
Jn 2 3 4 


Substitution for the integrals and using Equation 117 for the derivative of bulk 
axial velocity gives the final result 


K*,r)= CwV °M l 

r 1 */ + *v 


x44 + 4)4+44K 


RAs 4 


*5 

A 


(119) 


The parameters b, c, and d are known from the axial velocity solution of 
Appendix B so that Equation 1 19 can be readily solved. 


129 



Appendix D 

PRESSURE DETERMINATION 


The objective is to determine the axial pressure distribution in the liquid. 
The pressure distribution is required to' verify that the axial pressure drop does 
not exceed the maximum capillary pumping capability of the heat pipe. 

The axial momentum integral equation developed in Appendix B is 


±.\\ rdr+ ±. tv,*. 

Ot in CX in 

i\ v Au 


R f ~ R v \ dP 

dx 


+ v 


R, 


qu 

dr 


R, 


du 


R, v dr 


A 


+ C w R v ■ 


dV 


^ (120) 


0 


dx 


where the substitution has been made using the result derived in Appendix C 


d r R ' 

Wo = 17i R urdr 


Substitute the bulk axial velocity 


U = 


o r R t 
= urdr 

rJ - r; j r v 


and rearrange to give 
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1AL = _ AL 

^ dt dx 


2p 


rJ - R; ) dx J R 


2 

u n 

J D 


rdr + 


2 n 


1 J 

RT- r: 


R> 


du 

dr 


- R 


du 


R, v dr U, 


+ C w R v ■ 


avo 

dx 


Transforming to w* and r\ substitution of u* = br * + cr * 2 + Jr 


rearranging gives 


at/ _ _ dp 


at 


uf d^U 

“7 — f ^ 7 


+ >- u 


where 


y = 




r}-r;/ l 


/?/ /?/ 

a + 2c— +3</— . 

<5 5 J \ {Ri + R v ) 2 


4 


st? v 


K,v ( 


1 Tt( 3 + 4> 2 + 4 4 + 44 + M 5 + 44 


*v 


1 


R„ 


1 


R 


V \ .2 


+ -7- 5 + 6—- )c * + — 6+ 7-22 W+ 7 + 8— Jrf 
30 l 5 21 V 6 56 \ 5 



Use the trapezoidal rule to integrate Equation 122 


p[^+D_ £*")] = 


2 


dP 

dx 


+ fi- 


a 2 u 7T 

— r- + 

a* 2 


:(«+!) 


+ 


dP 

dx 


+ P 


d 2 U jj 
— r- + yU 

dx 2 


(n) 


Assume £/ ( ' 1+1) is known from continuity. Center difference the 
operator on U, 


d 2 u Ui-i-2V i+ U i+l 


dx 2 


(121) 


3 , and 


(122) 


(123] 

spatial 
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and rearrange to give 


dP ( n+ 1 ) 
dx 


_ 2 u ^\y~ k 


2 P 2fi \- , fi 


U i + 


~ U i+\ 


(«+l) 


+ 


(124) 


P 77 


u i-\ + \y + ~r~ 


2 P 2\x \- u - 


U ; +^rU i+l - 


dP 

dx 


|(») 


Difference the spatial operator on pressure. 


dP P i + 1 ~ P i 
dx 


to give 


p i+ =-[^. -pi* + 

n rr V . 2 P h x 2^\_ n - 


(n+1) 


+ 


(«) 


(125) 


All terms on the right hand side of Equation 125 are known so that the 
solution can proceed by specifying pressure at one location in the heat pipe. 
In general there is a location where the liquid and vapor have equal pressures. 
This point of matching pressure depends on the operating conditions. At 
moderate heat fluxes and high vapor pressure this point is at the end of the 
condenser. In this case, the maximum liquid pressure also occurs at this point. 
As the axial heat flux, and the vapor speed, increase, the point of matching 
pressure moves toward the beginning of the condenser. Prediction of the point 
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of matching pressure requires detailed information of the vapor flow dynamics, 
which is beyond the scope of this study. The point of matching pressure will 
be assumed to occur at the end of the condenser. With liquid pressure known 
at x = L, i = N, Equation 125 can be evaluated as 


p(rt+l)_ p {n+ 1) p{n) p ( 
r i - r i+ 1 + r /+ 1 - ‘i 


in) 


H — f 2 ph x 

+ {’>*- -IT 

u — f 2 ph x 


2 i n <" +1) 

h x) h x (126) 

2fi \— n — M 
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Appendix E 

ENERGY EQUATION NUMERICAL SOLUTION 


E.l GENERAL NUMERICAL FORMULATION 

The two dimensional energy equation in cylindrical coordinates neglecting 
viscous dissipation is 


dT , dT , cT 

— h u— b v— — = a 

at ax or 


d~T d^T 

~ 2 - 2 
C x or 


+ -LiI 

r dr 


The energy equation is to be solved subject to the boundary conditions 


AXIAL RANGE 

RADIAL RANGE 

CONDITION 

O 

II 

H 

R v < r < R w 

o 

II 

1 

x= 1 

Ry< r< R w 

O 

II 

^1 H 
1 

0 < x < L 

> 

II 

i 

Cu 

II 

0 < x < L e 

r=R w 

. dT _ 

*7 in 

L e < x < L e + L a 

II 

1*. 

k ar =o 

K dr ° 

L e + L a < x < L 

r=R* 

-kjl—mdl- 7J) 
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where subscripts u\ /, and v refer to the wall, liquid, and vapor respectively, 
q ir is the specified input heat flux, T e is the radiation environment rejection 
temperature, and q 0 is the liquid-vapor phase change heat flux given by 
Equation 30. 

Solution of the energy equation is accomplished by the finite difference 
Peaceman-Rachford Alternating Direction Implicit (ADI) method. The 
energy equation is rearranged as 


cT_ 

dr 



+ a 


d 2 T 
dx 2 



(128) 


Trapezoidal time integration is applied to Equation 128 to give in operator 
form 



+ aD L 

+ aD 0x 


+ 

+ 




(129] 


Like time levels are combined to give the Crank-Nicolson formulation 

{' - y[“A).r + + (f - v)d 0 , - uD 0 J} 7< W 11 = 

{ / + + “At + (t - v)c 0r - uD 0 J }^" 1 


The Crank-Nicolson formulation is second order accurate in space and time. 
The left hand side of Equation 130 can be factored in the form 
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/ - j-(«Dl r + (-2. - .-)£>o, r )][/ - y («£>L - uD 0jc ) 

I - y[<*A).r + aD 0 ,x + (f~ '’) D 0 .r - uD 0 J + 

'J 

A“( a£> 0 ,r + ( 7 " “ v )^0,r)( a ^0,.x _ uD 0,x) 


Similarly the right hand side can be factored as 


+ -r( aD 0 ,r + (7- “ V ) Z) 0,r)][ / + y( aZ) 0,x - uD 0 ,x) 
1 + ^"[ a ^0,r + a ^0^c + (7" ~ v ) £> 0,r“ w£) 0 ,x] + 

'y 

~r{ aD 0 \r + (l r ~ y ) D 0 ,r)( aD 0 ,x~ uD 0 ^) 


Substitution of the factored expressions in Equation 130 gives 


-K 


2 K aD 0s - uD ^ 


1 + “ uD 0^c 


' - y( aZ) 0,r + (t “ v ) £> 0,r) 

' + y( aZ) 0,r + (t"~ V ) D 0 ,r) 

J y{[( aZ) 0 ( r+ (t ~ V ) Z> 0,r)( aZ) 0^c ~ mZ) (u)] r ” + * “ 

[(«< + (t - »K)HL, - -Au)]^} 


jin+\) 


y 

)V n) 


+ 


(131) 


Equation 131 is consistent with the original Crank-Nicolson form to the last 
term on the right hand side. This term is the error of the factorization step. 
The factored equation can be shown in general to preserve the second order 
accuracy of the original Crank-Nicolson formulation. The 

Peaceman-Rachford time split of factored Equation 131 results in two coupled 
equations 
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/ - y («t>0,r + (t - '’Hr)] 7 " = [/ + j{ aD L - <‘t>0. x ) 


jin) 


1 ~ ~ uD °^) 


ji n + 1) _ 


1 + T ( aZ) 0 ,r + ( 7 - - v ) D 0,r) 


T 


where r* values are evaluated at an intermediate time t n < t* < t n+] . 
case the intermediate time is taken as r* = t n+Vl . 

Second order accurate centered differencing is used for the 
operators. 


^Ir = 


v»-iyy l +<>( ^ 


D o, x = 


T ,-\j- 2 T ij+T i+u 2 

0 + 0 (h x ) 


D 0,r = 


DjJ + 1 V l 

2K 


+ O(^) 




2/2. 


Substitution of the spatial operators, rearranging, and defining 


produces the final result 



( 132 ] 

( 133 ] 
In this 

spatial 


137 



where 


A 1 T'j_ | + A 2 T'j + A 3 T‘j + 1 - 

JK? + Aftr" + A «Tftj - + fi 5^ + VyH 




.4 2 = 0 r + 2a 




s, = 


a + w- 


J ^JC 


6 r 

B 2 = 6 r - 2a— 

^.X 


£i = 


a — u- 


°r_ 

j ^x 


A d = - a - u- 


A s = d x +2ct 


A 6 =-a + u-j 


B a = 




0 , 


J "r 


(134) 

(135) 
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Bs=0 x -2* 


6 r 


*6 = 




Equations 134 and 135 are sparse, tridiagonal matrices suitable for efficient 
solution by upper-lower decomposition and Gaussian elimination. The 
procedure is implemented by operating on the entire mesh with Equation 134 
to produce time t* temperatures. Equation 135 then sweeps the mesh of time 
t* values to produce time values. The only remaining issue is 

implementation of boundary conditions. 


E.2 IMPLEMENTATION OF BOUNDARY CONDITIONS 

E.2.1 Evaporator Endcap Boundary Condition 

x = 0 R v < r < R w 

The boundary condition is 

dx 

Since the solution derivative but not the solution is prescribed on the 
boundary, difference Equation 127 holds on the boundary. The boundary 
image point method is used. 

Center difference the boundary condition at the boundary i = 1 


139 



dT T oj 

dx 2 h x 

Solving for the image point T 0J gives 


Apply Equation 134 to the boundary 

A i T\j_ i + A 2 T*j + ^ 3 r^ +1 = Si7o^+ B 2 7 4 {j + B 2 7i n ] 
Substitute for the image point to give the first time split equation, 

AJ\j- 1 + A 2 r iJ+ A i r iJ+ 1 = * 2 7 #+ + * 3 ) t # 

In a similar manner, the boundary form of Equation 135 is 

A,T 17” + (A a + = B A r XM + fi 5 rr,. + B t r lJ+x 


E.2.2 Condenser Endcap Boundary Condition 

x-L R v < r < R w 


The boundary condition is 


dT 

dx 


= 0 


Center difference the boundary condition at the boundary i = N 


dT T N+l J T N-lj 

17~ 2 h x 


(136) 


(137) 
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The image point is T x+iJ , so that 

r v+l J= t N-\ j 


Following the same procedure as for the evaporator end cap boundary 
produces the two time split equations for the condenser endcap 

+ A ^j+ ^>1 = < B , + fi 3 )7 £> „ + B 2 7 $ 


K|+ -<6>b£lJ + A 5^ ]) =B A Tlj^ + B 5 r KJ +B 6 nj, 


l /+ 1 


E.2.3 Liquid- Vapor Interface Boundary Condition 

0 < x < L r = R., 

V 

The boundary condition is 


, ST, 

~ h— = 90 


where q 0 is given in Appendix A as 

The phase change velocity V 0 is determined in Appendix C from kinetic 
theorv as 


v o = 


j M 


/ 2 nR u 

Pi 

L 


- t s ) 


so that 
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^--c^[,4]T^-r s) 


or 


?0=Vv' 3,2 ( r v-r s ) 


The boundary condition becomes 


- */4^ = #X 3,2 (r v - t s ) 


or 


Center differencing at the boundary j = JV and solving for the image point 
Tijv-x gi ves 

t UV -> = t uv+ , + i^-r; 3/2 (r v - r yK ) 

Application of Equation 134 to the boundary j = JV and substituting for 'the 
image point gives 


2h r h 


a 2 ~ A \ 


r_q -.*—3/2 


T*JV + (^1 + A ?)^iJV + 1 ~ 


5 1 T*;”\ n/+ BoT^ n )i/+ B-xT^+x jy- A 


2h r h 


(138) 


i \ 1 i-\JV + *2 l \JV + ^1— 


4 v 


Note that in Equation 138, T* v is the vapor temperature at the intermediate 
time t * . This temperature is unknown a priori so that iterations are required 
between the energy equation solution, the flow dynamic equations, and the 
vapor determination. Since the unknown vapor temperature terms in 


142 


Equation 138 are nonlinear, special treatment is required. Quasilinearization 
is used to reduce Equation 138 to a relation with only linear unknown terms. 
The first nonlinear term is quasilinearized by 

( T* m+] )~ 3/2 - (7';™)- 3/2 + d( T* m y 3/2 [ T* m+ 1 - T* m ] 


or 

^ j.*m+l\-3/2 ^ 5 ^*my-3/2 3_^ 1 

d- 2 

where m is the iteration counter index. The above is now a linear expression 
in terms of T^ m+] . Similarly, the second nonlinear term becomes 

{r; m ^Y h & - ±-( T ; m )- 3l2 T ; m+ 1 

Substitution for the nonlinear terms in Equation 138 gives the final result for 
the first time split 


A'-) — A 


2h r h q r 


(tt)- 3 ' 2 - |( 7 rr 5/2 7 -; ra+l ]} 7 -ur+ m, + ^ K+1 


2 Kh, 


^1^-1 jv+ b 2A%+ 7$\ jv - A X — r -±\ Ur; m )-' A - 1 -[r; m )-^ 2 T; m+x 


/ L 


The second time split equation is found by applying Equation 135 to the 
boundary j= JV and substituting for the image point T iJV _ x 


JV + + A (A+\, JV = 


2 hh ‘ 

Bs ~ B 4 —^-T ;- il2 

K l 


2h r h. 1 .. 

T‘jv +(B a + B 6 )T*j V+1 + 71 

k l 
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Nonlinear terms in this second time split equation are for time t * values, which 
are known at time t {n+]) for the second split. Special treatment is not required 
for these terms. 


E.2.4 Pipe External Wall Boundary Conditions 

0 <x < L r — R w 


The boundary condition is 

i 5T . 

— ~ tfin Qrad 

where the input heat flux q in and the radiation heat flux q rad are functions of 
axial position 

Evaporator 0 < x < L e q in —* specified q rad = 0 

Adiabat L e < x < L e + L a qi n = 0 q rad = 0 

Condenser L e + L a < x < L q in = 0 q rad = — T\) 


Center differencing the boundary condition and solving for the image point 
T uw+ , gives 

T iJW+ 1 = T iJW- 1 “ \_Qin + zo {Tt ~ ^)] 

K w 
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where the product eo is taken as zero except for the condenser region. 
Application of Equation 134 to j = J\V and substituting for the image point 
gives for the first time split 


(A i + A 3 )T*j lv _ l + 


A 1 - A 


2h r e ° ^*3 
— T J UW 

2h r 


KjW ~ 


B l 7 i-\ ,JW + B 2 t]jw+ B i T \^\ y JW + Ay-j^Ji'iri - £cr T e A _ 


The first time split contains the unknown nonlinear term in wall temperature 
T ir nv ( the environmental radiation rejection temperature T e is assumed known 
and specified ). Quasilinearization is required 


(j*m+ 1)4^ \ 


Substitution in the first time split and rearrangement gives 
(A i + a 3 )T*j W _ 1 + 


8 h r e a 


4-4 r 

h 2 , 

\ 1 iJW) 


r r *m+ I 

liJW 


B \ J \-\jw + b iAjw+ B 3 T i+\jiv+ ^3-p L { < 7/« - Ea \}(T*jlv) A + T^ 4 ]} 

Application of the second time split equation to the boundary j = JW and 
substitution for the image point gives 


a aA^\JW + A 5 7 ) + A 6 = 

{ B 4 + B () r *iJW- \ + B S T hw~ B (q^~\ 4 in + £ff [ T ijw - T'e 4 ]} 


The nonlinear terms in the second time split are known since they are from 
time level r* . Special treatment is not required. 
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E.2.5 Corner Boundary Conditions 

Corner boundary conditions are evaluated by substituting for the image 
points across the boundary in each of the two boundary directions at the 
corner. The results are given below. 


EVAPORATOR END/EXTERNAL WALL CORNER 
The first time split equation is 


Ml + + 


%h r eo , m x 3 
A 2~ A 3 — Z 1 T \ JW ) 


T 


>m + 1 
\JW 


B 2^\ JW Jr Ml + ^ 3 )^ 2 n jW+ £<7 [ 3 (^iJh') + T e ]j 


The second time split equation is 

a 5 rfjiP + M4 + A 6^2 JW = 

M 4 + B 6 )T\jw- 1 + B 5 T \JW~ B 6~~r~\ q *in + ta \j\JW~ T e ]} 


CONDENSER END/EXTERNAL WALL CORNER 
The first time split equation is 


Mi + a 2) t NJW- 1 + 


Mi + b $I*n-\jw + 



T -*m + 1 
1 NJW ~ 



-*m 


NJW) 


X4 + r ; 4 


]} 
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The second time split equation is 


(A 4 + ^6) 7 A'-MH /+ A S^SJW ~ 

(£4 + B 6 )T* x jw _ j + B s T* N j W - B 6 -—^q* n + - T * 4 j 


EVAPORATOR END/LIQUID- VAPOR INTERFACE CORNER 
The first time split equation is 


a 2 ~ A r 


2h r h qV 


^ tn\ — 3/2 

kj L2^ T * ' 


1 + A \ 


- \{Tv m Y il2 T' m+ ' ] ^T\j V + (/I, + ^)7-,> +1 

-^p[y(7-;T' / ' - Yfn'T 3 ' 2 ^^'’ 


The second time split equation is 


A srfjv * + ( A 4 + ^6)^2JT^ = 


R 5 - B 4 - 


2A A T -m 


2h h 

T \JV + ( B 4 + + £ 4 — ]~~ T * 


Vz 


CONDENSER END/LIQUID- VAPOR INTERFACE CORNER 
The first time split equation is 


*1 - ^~[j<Trr il2 -}{rrr 5l2 T;^']y N j V+ (^ + 

(s, + b^_ uv+ B j&y- A ^.[±{rrr' A - |(rr)- 3 ' 2 7r +f 
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The second time split equation is 



T*N,JV + 1 + 


Appendix F 

THERMOPHYSICAL PROPERTY RELATIONS FOR SODIUM 


Selected thermophysical property relations are included in the program to 
update property values based on temperature conditions. All properties with 
the exception of surface tension are from Fink and Leibowitz[41 ] . The 
relation for surface tension variation with temperature is taken from 
Foust[42] . The relations that follow have been converted to consistent SI 
units. Property symbols and units are: 


Vapor Temperature 

T 

K 

Liquid Density 

Pi 

kglm 3 

Vapor Density 

Pv 

kgjm 3 

Liquid Thermal Conductivity 

k, 

Wfm K 

Thermal Diffusivity 

«/ 

m 2 js 

Specific Heat 

C P 

JlkgK 

Surface Tension 

a 

Njm 

Saturation Pressure 

P sal 

Pa 

Latent Heat of Vaporization 

kfg 

m g 

Dynamic Viscosity 

P 

kgjm s 
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The property relations also make use of the following identities: 

Critical Temperature T c = 2509.46 K 

Critical Density p c = 214.1 kgjm 3 

Reduced Temperature T r =T\T c 

Liquid Thermal Conductivity 

k { = 105.78 - 5.1767 x 10 -2 7"+ 4.8696 x 10“ 6 T 2 

Surface Tension 


a = 0.23402 — 1 x 10 V 


Saturated Liquid Vapor Pressure 


P ' . = 101,325 exp 


18.832- 13113 —1.0948 In T + 1.9777 x 10 V 


Latent Heat Of Vaporization 


h fg = 1.4487 x 10 6 (1 - T r )+ 3.4860 x 10 6 (1 - T r f 2 
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Dynamic Viscosity 


M = I.I259 x 10' 5 (p,) 1/3 exp 


0.74908 


Pi 


] 


Liquid Density 


Liquid density is dependent on the temperature ranges 


371 K < T< 1644 AT 


P/ = 101 1.8 - 0.22054 T — 1.9226 x 10 _5 7' 2 + 5.6371 x 10" y r 


9 - 7-3 


1 644 K < T< 2509 A: 


P/= P c [ 1 + 2.3709(1 - r r ) 0 - 31645 + 2.8467 x 10~ \t c - T ) 2 ] 


Vapor Density 


Pv = 


n fg 


1%) 


Pi 


sat 


-1 


which is the Clausius-Clapeyron Equation solved for vapor density with 
dP 


dT 


given under specific heat. 
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Specific Heat 



exp 18.832 - 13 ^ - -1.0948 In T+ 1.9777 x 10 4 T 


f dp A = 1.77587 x 10~ 10 

\ dp JT L b+cT r +dT 2 J 

- [2.15463 x 10 _21 r+ 1.220408 x \0~ 24 T 2 ]P sat + 1.671245 x \0~ 29 TP 2 

b = 1.2773035 
c = — 1.8267670 
d = 0.54946350 
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7< 1644 A' 



sat 


= 1.53139 x 10 3 - 0.6134527+ 3.35513 x 10 _4 7 2 + 


5.40593 x 10 6 



sal 


= -0.22054 - 3.8452 x 


io 5 r+ i.69i i x io _8 r 2 


7> 1644 A 


= 805.80 + 304.21(1 - 7 ) _ 0-67773 

01 sat 


= - -^-(.75027X1 - 7 r ) 0-68355 - ^5.6934 x 10~ 7 X7 C - T) 
0 1 J sat 1 c 


Thermal Diffusivity 


«/ = 


k l 

-Pl c P 
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Appendix G 

COMPUTER CODE LISTING 
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COODOOOG 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

COOOOOOO 

c 


S3CCQCOOCID30COOOOOOOCCDOC3CC3CCOOCC30COC00033C330C3CO: 

GREG ROCHE 
UCLA 

MECHANICAL. AEROSPACE AND 
NUCLEAR ENGINEERING 

PROGRAM HEAT. PIPE 

ANALYTICAL DETERMI NATI ON OF HEAT PIPE LIOUID DYNAMICS 

02 / 01/88 

00 0000000000 000 ooooooooooooooooooooooooooooooooooooooo 


0300000000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0000000000 


C 

c 

c 


IMPLICIT REAL#8 (A-H.O-Z) 

PARAMETER CIE=3. 

1 IC=6 . 

2 IMAX=21 . 

3 JV=1 . 

« JL=4, 

5 JW= 12 . 

6 IJMAX=IMAX*JU) 

COMMON /FLOW/UC I MAX, JV: JU) , VC IMAX, JV: JU) 

COMMON /FLOWS/US (IMAX. JV : JU) , VS ( IMAX, JV : JU) 

COMMON /FLOWN/UN (IMAX, JV : JU) , VN ( IMAX, JV : JU) 

COMMON /FLOWT/UT ( IMAX. JV : JU) , VT (IMAX, JV : JU) 

COMMON /UB AR/UB ( I MAX ) , GAMMA ( I MAX ) 

COMMON /TEMP 1/T ( I MAX , J V : JU ) , TN C IMAX, J V : JU ) , TUS ( I MAX ) 

COMMON /TEMP2/TV1 , TV2, TAV, TEA 

COMMON /TEMP3/TVES, TVESM, TVEP1 , TVCS, TVCSM, TVCP1 
COMMON /LUD/AC I JMAX, 3 ),BCI JMAX) 

COMMON / AD IASI/ ADI A1 ( JU) , ADI A2 ( JU) , ADI A3 ( JU) 

COMMON /ADIBS1/ ADIB1 ( JU) , ADIB2 ( JU) , ADIB3 ( JU) 

COMMON /AD I AS2/AD I A4 ( JU) , ADI AS ( JU) , ADIA6 ( JU) 

COMMON / ADIBS2/ ADIB4 ( JU) , ADIB5 ( JU) , ADIB6 ( JU) 

/ADIPAR/PHIR, PHIX, PHIXOR, PHIROX, PHIXRH, PHIRXH 
COMMON /F LUX/QS TR ( IMAX) , OIN 

COMMON /LPROP/RHOL, TKL, HFG, ALPHAL, SIGMA, UMU.UNU, CP 
COMMON /VPROP1/VMU , VNU , RHOV , PVAP 

SK rYIhoSS' RH0E ' RH0C - DPDTE • PS I E . fsie 

COfinON /UPROP/SRADC IMAX) , ALPHAU, TKU 

rnESnM /dx£t R n iYl : *2 C , J !^ J * RS1 ( JW 5 * RS2 C JU) * RS3 ( JU J * RS<i « JU 5 . RS5 ( JU) 

222222 ^5« PI ’ DEL TA * RV.RL.RU. JLfll . JLP1 , JVP1 . JVP2 , JVP3 
£222™ '2<p/x(inAX).iEm.iEPi.icm.icpi.inni.inn2 
connoN /xdim/el . al . cl . pi l . elph, clph 
common /deltas/hr. hx. ht. hxhalf, hrhalf. ktm 

COMHON /ITER/EPSTV.MAXITV. ITRTV 

connoN /uick/cu 

COMMON /PROPG/PG1 , PG2 . PG3 , PG4E. PGAC, PG5E . PG5C, PG6 , PG7 
COMMON /f=GROUP/F0 . FI . F2 . F3 . FA . F5 . F6 . F7 . F8. T13 
COMMON /CGROUP/C1 . C2 . C3 , C4 
COMMON /VGROUP/V1. V2, V3.VA.V5 

2?222o T /2GR0UP/21 . 22 . 23 . ZA . ZS . 26 . 27 . 28 . 29 . 21 0 . Z1 1 . 212 . 213 
DImInIiON ^^J* PLNCIr,AX) ' 2ETA(lr,AX) * 2ETAN finAX).RCAP(IMAX) 

DIMENSION TUC IMAX >,TS( IMAX) 

DATA ZETA, ZET AN/IMAXXO . DO . IMAXXO . DO/ 

DATA PL.PLN.RCAP/IMAXX0.D0, IMAXXO. DO. IMAXXO DO/ 

DATA TU.TS/IMAXX0. DO. IMAXXO. DO/ 

READ INPUT DATA 

READ (5.9010) NUM. (XPUT( I ) . 1 = 1 .NUM) 
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non non non non 


SET PROGRAM CONTROL CONSTANTS 

EPSDIT= 1 . 5D-1 
EPS0=5 . D-2 
EPSRG1= 1 . D-2 
EPSTS=1 . D-5 
EPSTV=1 . D-5 
EPSTW=1 . D-4 
MAXI TR=30 
MAXI TV=20 
ITFLTR=0 
IREGin=l 

NCASE=DINT CXPUT ( 1 ) ) 

SET AXIAL PARAMETERS 

HX=XPUT (2 ) 

HXHALF=HX#0.5D0 
IEM1=IE-1 
IEP1=IE+1 
ICM1=IC-1 
ICP1=IC+1 
IMM1=IMAX-1 
IMM2=IMAX-2 
PI L=HX*DFLOAT ( IMM1 ) 

EL=HX#DFL0ATCIEM1 ) 

AL=HX*DFLOAT< IC-IE) 

CL=HX*DFLOAT ( IMAX— I C) 

ELPH=HX* C DFLOAT (IE)— 0.5D0) 

CLPH=HXX (DFLOAT (IMAX— IE)— 0 .5D0) 

DO AO 1=1,1 MAX 

X( I )=HX*DFLOAT ( 1-1 ) 

CONTINUE 

SET RADIAL PARAMETERS 

PI =D AT AN Cl . DO )#4 . DO 
HR=XPUT (3 ) 

HRHALF=HR#0.5D0 
JVP1=JV+1 
JVP2=JV+2 
JVP3= JV+3 
JLM1=JL-1 
JLP1=JL+1 
RV=XPUT <4 ) 

RL=RV+HRXDFLOAT ( JL-JV) 

RW=RV+HRXDFLOAT ( JU- J V ) 

DELTA=RL— RV 
DO 50 J=JV,JUI 

R ( J ) =RV+DFLOAT ( J- J V ) *HR 
RS( J )=(R( J )-RV)/DELTA 

RSI ( J )=DELTA/C2 . DO XHX* ( DELTAKRS ( J >+RV ) ) 
RS2 ( J )=RS( J )*RS( J )/2 . DO 
RS3(J)=RS2CJ)*RSCJ)*2.D0/3.DO 
RS4 ( J )=RS3 ( J )XRS( J )X3 . DO/4 . DO 
RS5 C J )=RS4 ( J )*RS( J >#4 . DO/5 . DO 
50 CONTINUE 

DO 55 1=1. IMAX 
RCAPCI )=RV 
55 CONTINUE 

SET UICK PARAMETERS 

PR=XPUT ( 5 ) 

PD=2.D0XPR 
WR=XPUT ( 6 ) 

CW=PR*PR/C CPR+UR)#(PR+UR> > 
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c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


58 


C 

c 

c 


c 

c 

c 


PIPE FIXED VOLUPIES 

VOP=PI*RL*RL#PIL 
VOC=PI*RVXRV*PIL 
VOA= VOP-VOC 
VOLE=PI*RVXRV*ELPH 
VOLC=PIXRVXRVXCLPH 

SET TIDE PARAMETERS 

HT=XPUT ( 7 ) 

IF (NCASE.EQ.l) THEN 
TinE= 0 .D 0 
KTM=0 
ELSE 

READ (12) TIME.KTM 
ENDIF 
KTM0=KTM 
KSTEPS=A0000 
KOUIT=KTn+KSTEPS 

SET PRINT CONTROL PARAMETERS 

ISKIP1=100 

ISKIP2=2000 

IPRNT1=KTM+ISKIP1 

IPRNT2=KTM+ISKIP2 

IREPT1=1 

IREPT2=1 

IBEGIN=KTM0+1 

IEND=KOUIT-l 0 

SET RADIATION OPTICAL PARAMETERS 

EPS=0 . 7D0 
SIG=5. 67D-8 
EPSIG=EPS*SIG 
TE=XPUT (8) 

TE4=TE#*4 
DO 58 1=1,1 MAX 

IF CI.LT.IC) THEN 
SRADC I )=0 . DO 
ELSE 

SRADCI )=EPSIG 
ENDIF 
CONTINUE 

SET INPUT HEAT TRANSFER PARAMETERS 

0MAX=-1.5D6 
TAU=300 . DO 
OMEGA=PI/TAU 
JOFLTR=0 

01=2 . DOXPIKRWXHX 
02=2 . DOXPIKRVXHX 

Q3=PIXRVXRVXDSQRT(5 . D0*8 . 31AD3/(16 . D0X22 . 9898D0 ) ) 
Q4=PI*RVXRV/DSQRT (PD ) 

05=5 . 04D-8XPIXCRLXRL— RV*RV)/(PIL*PR) 

SET WALL PROPERTIES 

ALPHAW=5 . 3D-5 
TKW=1 . 10D2 
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c 

C SET CONSTANT SODIUM THERMOPHYSICAL PROPERTIES 

C 

TCR=250 9 . 46D0 
RH0CR=214 . IDO 
RG=8 . 3 14D3/22 . 9898D0 
VMU=1 .4D-5 

c 

C SET C FACTORS FOR MISCELLANEOUS CALCULATIONS 

C 

C1=CWKRVKHX/(RL*RL-RV#RV) 

C2=DSQRT ( 1 . DO/ (2 . D0*PI*RG) ) 

' C3=2.D0*HR*C2*CU 
C4=2.D0*HR/TKW 
C 

C SET V FACTORS FOR VAPOR TEMPERATURE CALCULATIONS 

C 

V1=2.D0XPI*RVXRV*PIXRV*RV 

V2=PI*RVXRV*DSQRT(5.D0*RG/16.D0) 

V3=HT*2.DO*PI*RV*CWXDSORT(1.DO/(2.DOXPI*RG))/VOLE 
V4=HT*2 . D0*PIXRVXCW*DSQRTC1.D0/(2.D0XPI*RG) )/VOLC 
V5=-HTXCW*RV*DSQRT(2.DO#PI/RG)/VOC 
C 

C SET F FACTORS FOR GRID VELOCITY CALCULATIONS (SUBROUTINE SPEED) 

C 

T13=0 . 2D0X ( 3 . DO— RV/CRL+RV) ) 

DELTA2=DELTA*DELTA 

DELTA3=DELTA*DELTA*DELTA 

FO=RV/DELTA 

Fl= CRL*RL-RV*RV )*3 . DO/RV 
F2=-( 0 . 5D0+F0 )*DELTA3 

F3=C6.D0+F0*12.D0)*0.4D0XDELTA2#CW*RV/(RL+RV)-3.D0*CRLXRL-RV*RV) 
F4=CRL*RL-RV*RV)*T13*0 . 5D0/RV-2. DOXCRL+RV) 

F5= (5 . D-2+FQ/12 . D0-T13*( 1 . DO/12. D0+F0/6 .DO ) )*DELTA3 
F6=-(12.DQ/lO5.DO+O.2D0XF0-T13*(O.2DO+O.4DOXF0)) 

1 *4 . D0XDELTA2*CWXRV/(RL+RV)-(RL*RL-RVXRV)XT13X0 . 5D0 

F7= CT13XT13* ( 1. D0/60 .D0+F0/30. DO )+(l. DO/168. D0+F0/1 05. DO) 

1 -T13*(2 - DO/105 . D0+F0/30 . DO ) )*4 . D0XDELTA2*CW*RV/C RL+RV ) 

F8=2 . DOXDELTA/CRL+RV) 

C 

C SET Z FACTORS FOR LIQUID PRESSURE CALCS (MAIN AND SUBR. SPEED) 

C 

ZO=RV/DELTA 

Z00=-8.D0*RV*CW/CCRL+RV)*CRL+RV)) 

Z1=2.DQ/CRLXRL-RVXRV) 

Z2=4.D0*RL/CDELTA*CRL*RL-RV*RV>) 

Z3=6.D0*RL/CDELTAXCRL*RL~RV*RV)) 

Z4=Z00*(3 . DO+4 . DOXZO )/12.D0 
Z5=Z00* (4 . DO+5 . DOXZO )/10 .DO 
Z6=Z0 OX (5 . DO + 6 . DOXZO )/15 . DO 
Z7-Z00X (5 . DO+6 . DOXZO )/30 .DO 
Z8-Z0 0XC 6 . DO + 7 . DOXZO )/21 .DO 
Z9=Z00X(7 . DO+8 . DOXZO )/56 .DO 
Z1 0=1 . DO/HX 
Zll=-2 . DOXHX/HT 
Z12=-2. DO/HX 
Z13=4 . DOXHX/HT 
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QUO 


c 

c 

c 


80 


85 


SET INITIAL SODIUM THERMOPHYSICAL PROPERTIES 


IF CNCASE.E0.1) THEN 
T0=XPUT C 9 ) 
TVEP1=T0 
TVCP1=T0 


CALL PROPS 
DUM=PVAP-SIGMA/RV 
DO 80 1=1,1 MAX 
PL (I >=DUM 
CONTINUE 
ELSE 

READ (12) 

READ (12) 

READ (12) 

READ (12) 

DO 85 1=1 


READ (12) 
CONTINUE 
ENDIF 


TVEP1 , TVCP1 , TAV 

S I GMA , TKL , HFG, RHOL , RHOV . UMU , UNU , VNU . CP . ALPHAL . PVAP 

PVE , PVC , RHOE , RHOC, DPDTE . DPDTC 

PG1 , PG2 , PG3 , PGAE. PGAC. PG5E. PG5C, PG6 . PG7 

IMAX 

PL ( I ) , RCAP( I ) , ZETAN(I ) 


SET ADI NUMERICAL PARAMETERS 


90 


C 

C 

c 


PHIR=2 . DOXHRXHR/HT 

PHIX=2 . DOXHXXHX/HT 

PHIXOR=PHIX/PHIR 

PHIROX=PHIR/PHIX 

PHIXRH=PHIXORXHRHALF 

PHIRXH=PHIROXXHXHALF 

DO 90 J=JV, JW 

IF (J.LT.JL) THEN 

ADIA1 (J)=-ALPHAL+ALPHAL#HRHALF/R(J) 

ADIA2 ( J )=PHIR+2 . DOXALPHAL 

ADI A3 ( J )=-ALPHAL-ALPHAL*HRHALF/R ( J ) 

ADIAA ( J )=-ALPHAL 

ADI A5 ( J )=PHIX+2 . DOXALPHAL 

ADIA6 ( J )=-ALPHAL 

ADI B 1 ( J )=ALPHAL*PHIROX 

ADIB2(J)=PHIR-2.DOXALPHALXPHIROX 

ADIB3 ( J )=ALPHAL*PHIROX 


ADIBA ( J )s 
ADIB5( J )= 
ADIB6 ( J )> 
ELSE 

ADIAKJ) 
ADIA2( J) 
AD I A3 (J ) 
ADIAA(J) 


(ALPHAL-ALPHALXHRHALF/RC J) )XPHIXOR 
PHIX-2 . DOXALPHALXPHIXOR 
(ALPHAL+ALPHALXHRHALF/R(J))xPHIXOR 

-ALPHAW+ALPHAWXHRHALF/R( J ) 

PHIR+2 . DOXALPHAW 
-ALPHAW-ALPHAWXHRHALF/R ( J ) 

-ALPHAU 


ADI A5 ( J ) =PHIX+2 . DOXALPHAW 

ADI A6 ( J )=-ALPHAW 

ADIB1 ( J)=ALPHAWXPHIROX 

ADIB2 ( J )=PHIR-2 . DOXALPHAWXPHIROX 

ADIB3 ( J )=ALPHAWXPHIROX 

ADI BA ( J )= ( ALPHAU- ALPHAWXHRHALF/R( J ) )XPHIXOR 
ADI B5 ( J )=PHIX-2 . DOXALPHAWXPHIXOR 
AD I B 6 ( J )= ( ALPHAW+ALPHAUXHRHALF/R ( J ) >XPHIXOR 
ENDIF 
CONTINUE 


DETERMINE INITIAL LIQUID AND VAPOR INVENTORIES 

IF (NCASE.EQ. 1) THEN 
WL=RHOLXVOA 
UV=RHOVXVOC 
WT=UL+WV 
ELSE 

READ (12) UL.UVjUT 
ENDIF 
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c 

c 

c 


c 

c 

c 


100 

110 


130 

140 


c 

c 

c 


c 

c 

c 


200 

210 


220 

230 


240 

250 


EVALUATE INITIAL OPERATING LIMITS 

PCMAX=2 . DO KS I GMA/PR 
QS0N=Q3*RH0E#HFG*DSQRTCTVEP1 ) 
QENT=Q4*HFGXDSQRT ( SIGMAKRHOV ) 
QCAP=Q5*SIGMA*HFG/UNU 

SET INITIAL VELOCITIES AND TEMPERATURES 

IF (NCASE.EQ.l) THEN 
DO 110 1 = 1 , IP1AX 
DO 100 J=JV,JW 
UCI, J)=0.D0 
VCI, J)=0.D0 
TCI, J) = T0 
CONTINUE 
CONTINUE 
ELSE 

DO 140 1 = 1 , I MAX 
DO 130 J=JV,JW 

READ (12) UCI. J). VCI, J). TCI, J) 
CONTINUE 
CONTINUE 
ENDIF 


READ ADDITIONAL STARTUP DATA 


IF CNCASE.NE.l) THEN 
READ (12) DELPL 

READ (12) QIN, QOT , QEV, OCO, QEVAP, QAX 
READ (12) ITER 
READ (12) IREGIM 
ENDIF 


PRINT SET-UP DATA AND INITIAL CONDITIONS 


WRITE (6,9003) 
WRITE (6.9004) 
WRITE (6,9002) 
WRITE (6,9015) 
WRITE (6,9013) 
WRITE (6,9014) 
WRITE (6,9031) 
WRITE (6.9007) 
WRITE (6,9005) 
WRITE (6,9030) 
WRITE (6,9006) 
WRITE (6,9011) 
WRITE (6,9027) 
WRITE (6,9017) 
DO 210 1=1. IE 


PIL , RW, RL, RV, VOP, VOC, VOA 
WL.WV.WT 

HR, HX, HT , PHIR, PHIX 


TIME, KTM 
ITER 
IREGIM 
TVEP1, TVCP1 
DELPL, PCMAX 

OI N . QOT , OEV . QCO . QEVAP . QAX , QSON . QENT , QCAP 
POOL 

S I GMA » TK L , HFG . RHOL . RHOV , UMU , UNU , CP , ALPHAL 


WRITE (6,9021) I, JV.UCI, JV), VCI, JV),TCI, JV),PL(I),RCAP(I) 
DO 200 J=JVP1,JW,4 


WRITE (6, 9019) I , J , UC I . J ) . VCI , J ). TCI , J ) 
CONTINUE 
CONTINUE 


DO 230 I=IEP1 , ICM1 

WRITE (6.9021) I. JV,U(I. JV).V(I.JV),TCI,JV),PL(I),RCAP(I) 
DO 220 J= JVP1 , JW, 4 

WRITE (6.9019) I . J . U( I , J ) , VCI . J ) . TCI . J) 

CONTINUE 

CONTINUE 

DO 250 I=IC. IMAX.6 


WRITE (6,9021) I, JV,U(I, JV),V(I, JV),TCI, JV),PLCI),RCAPCI) 
DO 240 J=JVP1 , JW, 4 

WRITE (6.9019) I , J , UC I . J ) . VCI , J) , TCI , J ) 

CONTINUE 

CONTINUE 
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260 


C 

C 

c 


c 

c 

c 


WRITE (6,9015) 

WRITE (6,9032) 

DO 260 1=1 , inAX 

WRITE (6.9033) I.ZETAN(I) 

CONTINUE 
WRITE (6.9015) 

BEGIN TIDE STEPS 

TSTAR=TI nE+0 . 5D0#HT 

TIME=TinE+HT 

KTM=KTM+1 

ITER=0 

UPDATE PROPERTY DEPENDENT ADI PARAMETERS 


505 

C 

C 

C 


510 


511 


515 


517 


DO 505 J=JV 
ADIAKJ) 
ADIA2( J ) 
ADI A3 ( J) 
ADIAA(J) 
AD I A5 ( J ) 
ADIA6( J ) 
ADIBKJ) 
ADIB2( J ) 
ADIB3 ( J ) 
ADIBA(J) 
ADIB5 ( J ) 
ADIB6 ( J )■ 
CONTINUE 


. JLM1 

=-ALPHAL+ALPHALXHRHALF/R( J ) 

=PHIR+2.D0#ALPHAL 

=-ALPHAL-ALPHAL*HRHALF/RCJ) 

=-ALPHAL 

=PHIX+2 . D0XALPHAL 

=-ALPHAL 

=ALPHALXPHIROX 

=PHIR-2 . DOXALPHALXPHIROX 

=ALPHALXPHIROX 

=(ALPHAL-ALPHALXHRHALF/R( J) jxphixor 
=PHIX-2 . DOXALPHALXPHIXOR 
=(ALPHAL+ALPHALXHRHALF/R(J))XPHIXOR 


UPDATE INPUT HEAT TRANSFER 


IF (TSTAR.LE.TAU) THEN 

D0 = 51^*i-i^lE*^ •DO-DSINCTSTARXOMEGA+O . 5D0XPI ) ) 

OSTRC I )=Q0 
CONTINUE 
ENDIF 

IF (JQFLTR.EQ.l) THEN 
DO 511 I=IEP1.IC 
OSTRCI )=0.D0 
CONTINUE 


OIEnl=0.25D0X(QSTR(IE-2)+2.D0KQSTRCIEni )+OSTR(IE)) 
? 5 °2J! (OSTR(IEm)+2 * D0>( OS TR (IE)+OSTR(IEPl)) )) 
OST R ( IElil )=QIEm TR C IE >+2 * D0 * QSTR ( 1EP1 >+OSTR C IEP1+2 i 
OSTR(IE)=QIE 
OSTRC IEP1)=QIEP1 
SUMQIN=0.D0 
DO 515 1 = 1. IflAX 

• IF (I.EO.l.OR.OSTRCI + D.EQ.O.DO) THEN 
SUf1QIN=SUMQIN+QSTR( I )X0 .5D0 
ELSE 


) 


SUI1QIN=SUMQIN+OSTR(I ) 
ENDIF 
CONTINUE 
ELSE 


SUf1QIN=0 .5D0X(QSTRC1 >+QSTRCIE)) 
DO 517 1=2, I EMI 

SUf10IN=SUf1QIN+OSTR(I ) 
CONTINUE 
ENDIF 

QIN=Q1XSUM0IN 
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UPDATE OPERATING REGIME AND AXIAL VAPOR MASS TRANSFER 

IF CDABS(QEVAP).GE.0.98DO*DABSCOSON)) THEN 
IF CIREGIM.EQ. 1) THEN 
IPRNT2=KTM- 
ENDIF 
IREGIM=2 

VMD0T=V2*RH0E*DSQRT C TVEP1 ) 

ELSE 

VMDOT=0 . DO 

IF ( IREGin . EQ . 2 ) THEN 
IREGIM=3 
ENDIF 

IF (IREGin. EO. 3) THEN 

IF ( DABS C TVEP1-TVCP1 ) . LT . EPSRG1 ) THEN 
IREGIM=1 
IPRNT2=KTM 
ELSE 

VMDOT=DABSC OEVAP )/HFG 
ENDIF 
ENDIF 
ENDIF 

PSIE=VMDOT*HT/VOLE 
PSI C=VMDOT*HT/VOLC 

SAVE PREVIOUS TIME TEMPERATURES 

DO 550 1=1, I MAX 
DO 54 0 J=JV,JW 
TNCI, J)=T(I, J) 

CONTINUE 

CONTINUE 

SAVE PREVIOUS TIME VELOCITIES 

DO 600 1=1, I MAX 
DO 590 J=JV,JL 
UNCI, J)=UCI, J) 

VNC I , J )=V ( I , J ) 

CONTINUE 

CONTINUE 

INITIAL GUESS TSTAR VELOCITIES AND VAPOR AND WALL TEMPERATURES 

TVES=TVEP1 
TVESM=TVEP1 
TVCS=TVCP1 
TVCSM=TVCP1 
DO 650 1=1,1 MAX 
TWSC I ) = T Cl , JW) 

DO 640 J=JV,JL 
US C I , J )=U C I , J ) 

VSCI , J ) = V C I , J ) 

CONTINUE 

CONTINUE 

BEGIN ITERATIONS BETWEEN ENERGY AND FLOW 
ITER=ITER+1 

RESET TEMPERATURES TO INITIAL VALUES 


DO 1050 1=1, I MAX 
DO 1040 J=JV,JW 
TCI, J)=TNCI, J) 
CONTINUE 
CONTINUE 
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c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


1090 

1100 


1150 


1200 


1250 


C 

C 

C 


SOLVE FIRST ADI TIP1E SPLIT 

ISPLIT=1 
CALL SPLIT1 

USE LU-DECOnPOSITION TO SOLVE THE FIRST ALGEBRAIC SPLIT SYSTEM 
CALL LORI 

LOAD TIME TSTAR TEMPERATURE RESULTS 

DO 1100 1=1,1 MAX 
DO 1090 J=JV,JW 

TCI, J )=B( J+JWXC 1-1 ) ) 

CONTINUE 

CONTINUE 

APPLY SMOOTHING FILTER IF REQUESTED 

IF (ITFLTR.EQ.l) THEN 
DELTMX=0 . DO 
DO 1250 I=IE, IC 
DO 1150 J=JV,JW 

DELT=DABS( 1.D0-T (I-1,J)/T(I,J)) 

IF (DELT . GT . DELTMX) DELTMX=DELT 
CONTINUE 

IF (DELTMX. GT.EPSDLT) THEN 
DO 1200 J=JV,JW 

TEMPI =0 .25D0*(T(I-2,J)+2.D0XT(I-1,J)+T(I,J)) 

TEMP2=0 .25D0*(T(I-1,J)+2.D0XT(I,J)+T(I+1,J)) 
TEMP3=0.25D0* (TCI. J)+2.D0*T(I+1, J)+T(I+2. J)) 

T ( 1-1 , J )= TEMPI 
T( I , J )=TEMP2 
T( 1+1 , J )=TEMP3 
CONTINUE 
GOTO 1250 
ENDIF 
CONTINUE 
ENDIF 


UPDATE VAPOR TEMPERATURES 

TVESM=TVES 
TVCSM=TVCS 

IF (IREGIM.EQ.l) CALL REGIM1 (ISPLIT) 

If (IREGIM.EQ.2) CALL REGI M2 ( ISPLIT) 
TVES=TV1 
TVCS=TV2 

SAVE T STAR ITERATION WALL TEMPERATURES 

DO 1400 1=1,1 MAX 
TWSCI )=T (I , JW) 

1A00 CONTINUE 
C 

c UPDATE TIME TSTAR VELOCITIES 
C 

CALL SPEEDCISPLIT) 

DO 1500 1=1,1 MAX 
DO 1450 J=JV,JW 
US ( I , J ) =UT < I , J ) 

VS(I, J)=VT(I. J) 

1450 CONTINUE 
1500 CONTINUE 


C 

C 

C 
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INITIAL GUESS TIME N+l VELOCITIES AND TEMPERATURES 

IF (ITER.EQ.l) THEN 
TVEP1=TVES 
TVCP1=TVCS 
DO 2100 1=1 . I MAX 
TS ( I )=T ( I » JV ) 

TLK I ) = T ( I , JU) 

DO 2050 J=JV.JW 
UCI, J)=U5CI. J) 

V(I, J)=VS(I, J) 

050 CONTINUE 

100 CONTINUE 
ENDIF 

SOLVE SECOND ADI TIME SPLIT 

ISPLIT=2 
CALL SPLIT2 

USE LU-DECOMPOS I T I ON TO SOLVE THE SECOND ALGEBRAIC SPLIT SYSTEM 
CALL LORI 

LOAD TIME N+l TEMPERATURE RESULTS 

DO 2300 J=JV.JW 

DO 2250 1=1,1 MAX 

T( I , j )=B ( I+IMAXXC J-l ) ) 

2250 CONTINUE 
2300 CONTINUE 

APPLY SMOOTHING FILTER IF REQUESTED 

IF (ITFLTR.EQ.l) .THEN 
DELTMX=0 . DO 
DO 2500 I=IE, IC 
DO 2350 J=JV,JW 

DELT=DABS( 1 . D0-T (I-1,J)/T(I,J)) 

IF (DELT.GT.DELTMX) DELTt1X=DELT 
2350 ... CONTINUE 

IF ( DELTMX . GT . EPSDLT ) THEN 
DO 2400 J=JV,JW 

TEMP1= 0 . 25D0M ( T ( 1-2 . J )+2 . D0*T ( 1-1 , J >+T ( I , J ) ) 
TEMP2=0.25D0*CT(I-l.J>+2.DO#TCI,J)+TCI+l, J)) 
TEMP3=0.25D0#(T(I,J>+2.DOKTCI+l,J)+T(I+2, J)) 

TCI-1. J)=TEMP1 
T( I , J )=TEMP2 
TCI+1, J)=TEMP3 
2400 CONTINUE 

GOTO 2500 
ENDIF 

2500 CONTINUE 
ENDIF 

UPDATE VAPOR TEMPERATURES 

IF ( IREGIM. EQ. 1 ) CALL REGIM1 (ISPLIT) 

IF (IREGIM. EQ. 2) CALL REGIM2CISPLIT) 

TVEP1=TV1 

TVCP1=TV2 
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c 

C DETERMINE TIME N+l VELOCITIES 
C 

CALL SPEEDC ISPLIT) 

DO 2650 1=1,1 MAX 
DO 2600 J=JV,JW 
U(I, J)=UT(I, J) 

V(I, J)=VT(I,J) 

2600 CONTINUE 
2650 CONTINUE 
C 

C FIND ITERATION MAX NORM CHANGE OF TN+1 LIQUID SURFACE TEMPERATURE 

C 

DELTSM=0 . DO 
DO 2700 1=1, IMAX 

delts=dabs ( 1 . DO-TSC I )/T ( I , JV) ) 

IF (DELTS. GT.DELTSM) THEN 
DELTSM=DELTS 
I TS=I 
END IF 
2700 CONTINUE 
C 

C FIND ITERATION MAX NORM CHANGE IN PIPE WALL TEMPERATURE 

DELTWM=0 . DO 
DO 2750 1=1, IMAX 

DELTW=DABS ( 1 . DO-TWC I )/T(I , JW) ) 

IF (DELTW. GT.DELTWM) THEN 
DELTWM=DELTW 
ITW=I 
ENDIF 
2750 CONTINUE 
C 

C UPDATE TIME N+l ITERATION WALL AND SURFACE TEMPERATURES 


DO 2800 1=1, IMAX 
TS( I )=T ( I , JV) 

TWC I )=T ( I , JW) 

2800 CONTINUE 
C 

C CHECK FOR CONVERGENCE 
C 


C 

C 

C 

2850 

C 

C 

c 


IF ( I TER . GE . MAXI TR ) THEN 

IPRNT2=KTM° 26) ITER * EPSTS * DELTSf1 * ITS * EPST W.DELTWM, ITW, IREGIM 

KST0P=2 
GOTO 2850 
ENDIF 

iL,4 DELTSP1 - GT - EPSTS - 0R - DELTWri - GT - EPST wJ GOTO 1000 
CONTINUE 

UPDATE THERMOPHYSICAL PROPERTIES 
CALL PROPS 

UPDATE HEAT TRANSFER LIMITS 


PCMAX=2 . DOKSIGMA/PR 
QSON=Q3#RHOE*HFGKDSQRT ( TVEP1 ) 
QENT=Q4#HFG*DSQRT(SI6MA#RH0V ) 
Q C AP= Q5 * S I GMA * HFG/UNU 
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EVALUATE HEAT TRANSFER 

SUMOOT=0 .DO 
SUMQEV=Q . DO 

sunoco=o .do 

SUMQEI=0.D0 
DO 3000 1=1,1 MAX 

PHASE CHANGE HEAT TRANSFER 

QPC=PG3*V(I. JV) 

IF ( I . EQ . 1 . OR . I . EQ . I MAX ) QPC=QPC*0 . 5D0 
IF CVCI, JV).LT.0.D0) THEN 
SUMQEV=SUMQEV+0PC 
ELSE 

suriQCO=sunoco+OPC 

ENDIF 

IF (I.LE.IE) SUMOEI=SUMQEI+QPC 

OUTPUT HEAT TRANSFER 

IF (I.GE.IC) THEN 

QOI =EPS*SI G* (TCI, JW ) HKA-TE4 ) 

IF (I .EQ. IC.OR. I .EQ. IMAX) QOI=QOI#0 . 5D0 

sumqot=sumoot+ooi 

ENDIF 
000 CONTINUE 

OOT=Ql#SUriQOT 

QEV=Q2*SUMQEV 

QCO=Q2#SUMQCO 

QEVAP=Q2*SUnOEI 

QAX=VMDOT#HFG 

CHECK CAPILLARY LIMIT HEAT TRANSFER 

IF CDABS(OEV) .GT.QCAP) THEN 
WRITE (6,9000) QEV,QCAP 
IPRNT2=KTM 
KST0P=2 
ENDIF 

CHECK HEAT TRANSFER FOR STEADY STATE 

QTEST=DABS( 1 . DO+OOT/OIN) 

IF (QTEST.LT.EPSQ) THEN 
WRITE (6,9009) 

KSTOP=2 

IPRNT2=KTM 

ENDIF 

UPDATE LIQUID PRESSURES 

DO 3300 1=1, I MAX 
PLN ( I )=PL ( I ) 

3300 CONTINUE 

PL ( I MAX ) =PVAP-SIGMA/RV 
DO 3A00 I=IMM1,2,-1 

ZETA(I )=Z10*UMU#(UB( 1-1 HUBCI+l) ) 

1 + ( GAMMA ( I ) KHX+Z1 1 *RH0L+Z12#UMU >#UB ( I ) 

PL (I )=PL(I+1)+PLN(I+1)-PLN(I)-ZETA(I >-ZETAN(I> 
ZETAN( I )=ZETA( I )+Z13*RH0L*UB ( I ) 

3A00 CONTINUE 

ZETA ( 1 ) =- ( UB ( 3 ) -2 . D0#UB ( 2 ) ) KUMU/HX 
PLC1)=PL(2)+PLN(2)-PLN(1)+ZETA(1)+ZETANC1) 

ZETAN( 1 )=ZETA( 1 ) 
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FIND CAPILLARY RADII OF CURVATURE 


DO 3500 I=1,IMAX 

£^ A £ = C-3 - do *V(I»JV) + 4. D0*V( I , JVP1 )-V ( I , JVP2 ) )XUMU/HR 
350 0 CONTINUE 1 ^ ~ 2 ' D0 * SI GP1A/C PVAP_PL ( 1 >+PG7*V( I , JV )*V( I , J V)+GRAD ) 


CHECK CAPILLARY LIMIT 

DELPL = PL(iriAX)-PL(l) 

IF ( DELPL . GT . PCflAX) THEN 
. WRITE (6.9001) DELPL. PCMAX 
IPRNT2=KTM 
KST0P=2 
ENDIF 


HASS INVENTORY ACCOUNTING 

POOL=WT-RHOE*VOLE-RHOC*VOLC-RHOL*VOA 

PRINT NEW DATA 

IF ( KTfl . GE . IEND ) IPRNT2=KTM 
IF (KTn.EQ. IPRNT1) THEN 
IF (IREPTl.EO.l) THEN 
IPRNT1=IPRNT1+1 
IREPT1=0 
ELSE 

IPRNT1=IPRNT1+ISKIP1— 1 
IREPT1=1 
ENDIF 

IF (KTn.NE. IPRNT2) THEN 
WRITE (6,9022) 

WRITE (6.901A) TIME, KTM 
WRITE (6,9031) ITER 
WRITE (6,9007) IREGIM 
WRITE (6,9005) TVEP1.TVCP1 
WRITE (6,9030) DELPL, PCflAX 

cf ' 9011 ) POOL 001 ' ° EV * 0C ° ' 0EVAP ' 0AX ' 0S0N ' CENT , OCAP 

ENDIF 

ENDIF 

IF (KTM. EQ. IPRNT2 . OR. KTM. LE. I BEGIN) THEN 
IF (KTM.EQ. IPRNT2) THEN 
IF CIREPT2.EQ.1) THEN 
IPRNT2=IPRNT2+1 
IREPT2=0 
ELSE 

IPRNT2=IPRNT2+ISKIP2-1 

IREPT2*1 

ENDIF 

‘ENDIF 

WRITE (6,9015) 

WRITE (6.901A) TIME, KTM 
WRITE (6,9031) ITER 
WRITE (6,9007) IREGIM 
WRITE (6,9005) TVEP1,TVCP1 
WRITE (6,9030) DELPL, PCMAX 

WRItI (6*9011) pool OOT,OEV ' OC °' OEVAP ' OAX ' OSON ' OENT ' OCAP 
WRlH (6*9017) G ^* A * TKL , HFG, RHOL , RHOV , UMU , UNU , CP , ALPHAL 
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c 

C EVAPORATOR DATA 

C 

DO 5020 1=1, IE, 1 

WRITE (6,9021) I , JV, U ( I . JV ) , VC I . JV) , T ( I . JV ) , PL ( I ) , RCAP C I ) 
DO 5000 J=JVP1 , JL , 1 

WRITE (6,9019) I , J , U( I , J ) , V( I . J ) , T( I . J ) 

5000 CONTINUE 

DO 5010 J= J L+2 , JW, 2 

WRITE (6,9018) I,J,TCI,J) 

5010 CONTINUE 

5020 CONTINUE 

C 

C ADIABAT DATA 

C 

DO 5050 I=IEP1 , IC+3 , 1 

WRITE (6,9021) I , JV, UC I , JV) , VCI , JV) , T( I , JV ) , PL ( I ) , RCAP C I ) 
DO 5030 J=JVP1 , JL, 1 

WRITE (6,9019) I , J , UCI , J ) , VCI , J ) , T( I . J ) 

5030 CONTINUE 

DO 5040 J=JL+2, JW, 2 

WRITE (6,9018) I.J.TCI.J) 

5040 CONTINUE 

5050 CONTINUE 

C 

C CONDENSER DATA 

C 

DO 5080 I=IC+4, IMAX-3,3 

WRITE (6,9021) I , JV, UC I , JV) , VCI . JV) , T( I , JV ) . PL ( I ) , RCAPC I ) 
DO 5060 J= JVP1 , JL , 1 

WRITE C 6 , 9019 ) I , J . UC I , J ) , VC I . J) , TCI . J ) 

5060 CONTINUE 

DO 5070 J= JL+2 , JW, 2 

WRITE C6.9018) I,J,TCI.J> 

5070 CONTINUE 

5080 CONTINUE 

DO 5095 I = I MAX— 2 , I MAX 

WRITE C 6 , 9021 ) I , JV, UC I , JV) , VCI . JV) . TCI , JV) , PL C I ) , RCAPC I ) 
DO 5085 J=JVP1 , JL, 1 

WRITE C6, 9019) I , J, UC I . J ) , VCI , J ) . TCI , J ) 

5085 CONTINUE 

DO 5090 J=JL+2, JW, 2 

WRITE C6, 9018) I.J.TCI.J) 

5090 CONTINUE 

5095 CONTINUE 

WRITE C6.9015) 

ENDIF 

C 

C DO NEXT TIME STEP IF TIME HAS NOT EXPIRED 
C 

IF CKTM. LT . KQUIT. AND. KSTOP . NE . 2) GOTO 500 
WRITE C6.9008) TIME 
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c 

c 

c 


WRITE TO STARTUP FILE 


WRITE (13) TIME.KTM 
WRITE (13) TVEP1.TVCP1, TAV 

{jpy luf n m,Z K ii^ G 4, RH0L ' RH0V ' UPIU * UNU * VNU * CP » ALPHAL , P V AP 

r^iI E <13) pve.pvc,rhoe,rhoc,dpdte.dpdtc 

DO 7 02 0 1 1 i i ?I riAX G2 ' PG3 ' PGAE ' PGAC * PG5E * PG5C. PG6 . PG7 

7020 CONTINUE P '-CI ) .RCAPtI >.ZETAN(I ) 

WRITE (13) WL.WV.WT 
DO 70A0 1 = 1 , iriAX 
DO 7030 J=JV, JW 

7030 CONTINUE U3> 

70A0 CONTINUE 

WRITE (13) DELPL 

WRItI (13) f^OOT.OEV.OCO.OEVAP.OAX 

WRITE (13) IREGIP1 
8000 STOP 


FORMAT STATEMENTS 


C 

c 
c 

9000 FORMAT (/////22X, * XXKXX CAUTION XKKKK* 

•vaoZd L uc!t H to T transfer LIMIT EXCEEDED* 

2/16X, ’VAPOR NEAT TRANSFER = '.1PD12.5 

onn-i ^^Zn«i-!- CAP ^ 1 '*' ARY HEAT TRANSFER LIMIT = * , 1PD12 5) 

9001 FORMAT (/////22X, * XXXKX CAUTION XXKXX* 

iKZJS* CAPILLARY PUMPING LIMIT EXCEEDED* 

2/16X. * LIQUID PRESSURE DROP = * , 1PD12 5 
onno 3 cnn^* * MAX CAPILLARY PRESSURE = * , 1PD12 5) 

9002 FORMAT (/IX, * NUMERICAL PARAMETERS : * 

K^? X *’ RADIAL STEP SIZE (M) = ’ , 1PD10 .2 
2/13X, ’AXIAL STEP SIZE (M) = *,1PD10 2 
3//1AX. 'TIME STEP SIZE (S) = *!lPD102 
A//1X, ’ADI NUMERICAL PARAMETERS:* 

5/28X, * PHIR = ’ , 1PD10 . 2 
„„ 6/28X, ’ PHIX = * , 1PD1 0.2) 

9003 FORMAT C*l», IX, » h EAT PIPE DIMENSIONS:* 

1//13X, INTERNAL LENGTH (M) = *,D12 5 
2/16X, ’ OUTER RADIUS (M) = *,D12.5 
3/13X, * INTERNAL RADIUS (M) = *,D12 5 
A/17X, ’WICK RADIUS (M) = ',D12.5 

5/AX, 'TOTAL INTERNAL VOLUME (MXK3) = *,D12 5 
6/8X, VAPOR CORE VOLUME (MXX3) = '.D12.5 
7 ~£* ’ LIQUID ANNULUS VOLUME (MXX3) = *Id12 5) 

• Liooio T n!5» Z «s>* A I. m 2 E s“’" ,c1nb Am. masses:. 

2/17X, 'VAPOR MASS (KG) = *,D12 5* 
onnc 3 ^ii X * ’ T0TAL FLUID ^SS (KG) = * * ,D12.5) 

9005 FORMAT (1 OX, 'EVAPORATOR VAPOR TEMPERATURE (K) = ' 3PD20 14 

• o a * l ^:; C 9?25 NS f R VAP0R TE »P ERA TURE (K) = *.3PD20.l4f 
9006 s , :?oD A T <19X * 'INPUT HEAT TRANSFER (W) = * . 1PD12 5 
1/18X, ’OUTPUT HEAT TRANSFER (W) * * . 1PD12 5 

IxZoS' ! EVAP0RATI0N HEAT TRANSFER (W) = * . 1PD12 5 
9° NDENSATI0N HEAT TRANSFER (W) = * , 1PD12 5 
A/ AX, 'EVAPORATOR HEAT TRANSFER (W) = *,1PD12 5 
5/19X, ’AXIAL HEAT TRANSFER (W) = *,1PD12 5 ’ 

HEAT TRANSFER LIMIT CW) * *IlPD12.5 
7/7X ENTRAINMENT HEAT TRANSFER LIMIT CW) = *.1PD12 5 
on «-s 8 Z 9X * CAPILLARY HEAT TRANSFER LIMIT CW) = *,1PD12 5) 

F0RriAT C21X, ’OPERATING REGIME FLAG = ’.15) 12 

cS« MAT (/1X '’CASE COMPLETED' ,D12. 5) 

9010 SSKS ( 6X^1 2/ ( A ^6X?D12°6 ) ) ) "^ ** CASE C °" PLETED ****’> 

9013 FOR^J E/JiS: 'l^L^Dm^NS: V) E ** 1PD12 - 5) 
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9014 FORriAT (10X, 'xxxxx TIME = ' , 1PD12 . 5 . 5X, ' TIME STEP = ',18,' xxxxx*) 

9015 FORMAT (1H1) 

9017 FORMAT (/6X, * I ' , 7X, ' J ' , 12X, ' U’ , 16X, ' V ’ , 

+ 15X, ’ TEMPERATURE’ ) 

9018 FORMAT (AX, 13 , 5X, 13 , 42X. 1PD17 . 10 ) 

9019 FORMAT ( AX, 13 , 5X, 13, 2(6X, 1PD12 . 5) . 6X, 1PD17 . 10 ) 

9021 FORMAT ( AX, 13, 5X, 13 , 2( 6X, 1PD12 , 5 ) , 6X, 1PD17 . 10 , 6X, 1PD12 . 5, 
16X.1PD12.5) 

9022 FORMAT (///) 

9026 FORMAT (/////20X, ’ xxxxx CAUTION xxxxx* 

1/12X, ’ CONVERGENCE WAS NOT REACHED IN ITERATIONS = ',15 
2/15X, 'SURFACE TEMP CONVERGENCE TOLERANCE EPSILON = '.1PD12.5 
3/9X, ' MAX SURFACE TEMPERATURE NORMALIZED CHANGE = '.1PD12.5 
A/33X, ' INTERFACE GRID POINT = '.15 

5/12X, * WALL TEMPERATURE CONVERGENCE TOLERANCE EPSILON = '.1PD12.5 
6/15X, ' MAX WALL TEMPERATURE NORMALIZED CHANGE = ’.1PD12.5 
7/38X, 'WALL GRID POINT = '.15 
8/A2X, ' IREGIM FLAG = ’,15) 

9027 FORMAT (15X, 'XXXXX THERMOPHYSICAL PROPERTIES XXXXX* 

2/9X, 'SURFACE TENSION COEFFICIENT SIGMA = ’.1PD12.5 
3/1 IX, ’ LIQUID THERMAL CONDUCTIVITY TKL = ’.1PD12.5 
4/1 IX. ’ LATENT HEAT OF VAPORIZATION HFG = ’.1PD12.5 
5/23X, ’ LIQUID DENSITY RHOL = ’.1PD12.5 

6/24X, ' VAPOR DENSITY RHOV = '.1PD12.5 

7/14X, ' LIQUID DYNAMIC VISCOSITY UMU = '.1PD12.5 

8/12X, ’ LIQUID KINEMATIC VISCOSITY UNU = ’.1PD12.5 

9/19X. ' LIQUID SPECIFIC HEAT CP = '.1PD12.5 

1/9X. 'LIQUID THERMAL DIFFUSIVITY ALPHAL = ’.1PD12.5) 

9030 FORMAT (17X, 'LIQUID PRESSURE DROP (PA) = '.1PD12.5 
1/15X. 'MAX CAPILLARY PRESSURE (PA) = ’.1PD12.5) 

9031 FORMAT (10X, 'TIME STEP CONVERGENCE ITERATIONS = '.15) 

9032 FORMAT (/6X, ' I ' , 11X, 'ZETAN' ) 

9033 FORMAT (4X, I3.5X, 1PD12.5) 

END 
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x * BLOCK DATA INITIALIZATION * 

C ****£*J£***J*J**********#********#*KX*KKXKX«#*K####KK*KK###k##*kkk**## 

IMPLICIT REAL*8 (A-H.O-Z) 

PARAMETER (IE=3, 

1 IC=6, 

2 IMAX=21 , 

3 JV=1 , 

* J L=A , 

5 JW=12, 

6 I JMAX= I MAX* J W ) 

9- SD«2 N /F 1 01 U/U C I nAX * J V : J W ) , V C I MAX . J V : J W ) 

^LOUS/USCIMAX, JV : JW) , VS( IMAX, JV: JW) 

£ L OWN/UN ( I MAX , J V : JW ) , VN C I MAX , J V : J W ) 

/ ,f,9 0UJ/U1 ( IP1AX > JV : JW ) , VT ( IMAX, JV : JW) 

COMMON /UBAR/UB(IMAX),GAMMA(IMAX) 

COMMON /FLUX/QSTRCIMAX) , QIN 

DATA°U , V/I JMAXXO ^Dof I JMAXXO ! D0/ > ' RS2 ( JW) * RS3 < JW) , RSA ( JW) . RS5 ( JW) 
DATA US, VS/I JMAXXO . DO , I JMAXXO . DO/ 

DATA UN, VN/I JMAXXO . DO, I JMAXXO .DO/ 

DATA UT , VT/I JMAXXO . DO, I JMAXXO DO/ 

DATA UB/IMAXXO . DO/ 

DATA GAMMA/IMAXXO . DO/ 

P*I A DC7 S Dc^ 1 o!?f ^ 2/JWXO . DO , JWXO . DO, JWXO .DO, JWXO . DO/ 

DAT A RS3 , RSA , RS5/JWX0 . DO , JWXO . DO, JWXO DO/ 

DATA OSTR/I MAXX 0 . DO/ 

END 
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c XXXXXXXXKX#K#XXKXXXXXXXXXXXXXXX#XXXXXX#XXXXXXXXXXXXXXXX*K##K*#KXKK#KX# 

C X SUBROUTINE SPLIT1 - ENERGY EQUATION SOLUTION x 

C x x 

C X THIS ROUTINE USES PEACEP1AN-RACHFORD TEMPORAL ADI TO SOLVE THE X 

C X FIRST TIME SPLIT OF THE ENERGY EQUATION X 

C X * 

c * * 

c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

c 

c 


SUBROUTINE SPLIT1 
IMPLICIT REALK8 CA-H,0-Z) 
PARAMETER (IE=3, 


1 

2 

3 

A 

5 

6 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


I C=6 , 

IMAX=21 $ 

JV=1, 

JL=A, 

JW=12, 

I JMAX=IMAX*JW) 

/FLOW/UCIMAX, JV: JW), VCIMAX, JV:JW) 

/F LOUS/US C I MAX, JV : JW) , VS ( I MAX* JV : JW) 

/FLOWN/UN C I MAX, JV : JW) , VN( I MAX, JV: JW) 

/FLOWT/UTCIMAX, JV: JW), VTCIMAX, JV: JW) 

/UB AR/UB ( I MAX ) , GAMMA ( I MAX ) 

/TEMP1/T Cl MAX, JV : JW) , TN( IMAX, JV : JW) , TWSC IMAX) 

/TEMP2/TV1 , TV2 , TAV , TEA 

/TEMP3/TVES, TVESM, TVEP1 , TVCS, TVCSM, TVCP1 
/LUD/AC I JMAX, 3 ) , B ( UMAX) 

/ADI AS 1/ADI A1 C JW) , ADI A2( JW) , ADI A3C JW) 

/ADIBS1/ADIB1 C JW) , ADIB2C JW) , ADIB3C JW) 

/ADI AS2/ADI AA( JW) , ADI A5( JW) , ADIA6C JW) 

/ADI BS2/ADIBA ( JW) , ADIB5 C JW) , ADIB6 C JW) 

/ADIPAR/PHIR, PHIX, PHIXOR, PHIROX, PHIXRH, PHIRXH 
/FLUX/OSTR C I MAX ) , 0 I N 

/LPROP/RHOL , TKL , HFG, ALPHAL , SIGMA, UMU, UNU , CP 
/VPROP1/VMU , VNU, RHOV , PVAP 

/VPR0P2/PVE, PVC, RHOE, RHOC,DPDTE,DPDTC,PSIE,PSIC 
/CRI T/TCR, RHOCR 
/WPROP/SRADC IMAX) , ALPHAW, TKW 

/RAD/RC JW) , RSC JW) ,RS1 ( JW) ,RS2CJW) , RS3C JW) . RSA ( JW) , RS5 ( JW) 

/R/PI , DELTA, RV , RL , RW, JLM1 , JLP1 , JVP1 , JVP2 , JVP3 

/XD/XC IMAX) , IEM1 , IEP1 « ICM1 , ICP1 , IMM1 , IMM2 

/XDIM/EL , AL, CL# PIL, ELPH, CLPH 

/DE L T AS/HR , HX , HT , HXH ALF , HRHALF , KTM 

/I TER/EPSTV , MAXI TV , ITRTV 

/WICK/CW 

/PROPG/PG1., PG2 , PG3 , PGAE , PGAC, PG5E , PG5C, PG6 , PG7 
/FGROUP/FO , FI , F2 , F3 , FA , F5 , F6 , F7 , F8 , T13 
/CGROUP/C1 , C2 , C3 , CA 
/VGROUP/V1 , V2 , V3 , VA , V5 

/ZGROUP/Z1 , Z2, Z3, ZA,Z5,Z6,Z7,Z8,Z9,Z10,Z11 ,Z12 , Z13 


C 

C EVALUATE INTERMEDIATE TIME SPLIT MATRIX ELEMENTS 
C 


DO 199 1=1, IMAX 
DO 190 J=JV,JW 
K=J+JW*CI-1) 


C 

C EVALUATE COEFFICIENTS 

C 

A1=ADI A1 C J )-VS( I , J )#HRHALF 
A2=ADIA2 ( J ) 

A3=ADI A3 ( J )-VS( I , J )*HRHALF 
B1=ADIB1 ( J )+UN( I , J )*PHIRXH 
B2=ADIB2C J ) 

B3=ADIB3( J )-UN(I , J )#PHIRXH 


172 



c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


1 


c 

c 

c 

c 

c 

c 


c 

c 

c 


EVAPORATOR END BOUNDARY CONDITIONS 
IF CI.EO.l) THEN 

CORNER BETWEEN L/V INTERFACE AND EVAPORATOR END 
IF (J.EQ.JV) THEN 

TVTRm = 2.5D0/(TVESri#*l .5>-l . 5D0/ ( TVESn#*2 . 5 >*TVES 
TVTRH 2 - 1 . 5D0/DSGRT ( TVESID-0 . 5D0/( TVESn##l 5 )*TVES 
A(K,1)=0.D0 

AC K, 2 ) = A2-Al#PG<4EKTVTRni 
A ( K , 3 )=A1+A3 

B(K)=B 2 KT(I , J)+(B1+B3)#TCI, J >-Al*PG4E*TVTRn2 
ENDI F 

evaporator end between corners 

IF (J.NE. JV.AND. J.NE. JW) THEN 
A(K, 1 )=A1 
A( K , 2 )=A2 
A(K, 3 )=A3 

BCK)=B2*TC1, J)+(B1+B3)XT(2, J) 

ENDIF 


CORNER BETWEEN EVAPORATOR END AND AXIAL WALL 

IF (J.EQ.JW) THEN 
A(K, 1 )=A1+A3 

A( K, 2 >=A2-A3*4 . D0XC4XSRADCI )*CTWS( I )*#3 ) 

A ( K, 3 )=0 .DO 

BCK)=B2XT(I, J)+(B1+B3)XT(I+1, J) 

cunTC +A3XCA* ( OSTR ( I )-SRAD( I )* (3 . DO* ( TWS C I >#*<♦ ) + TE4 ) ) 
cNDIF 
ENDIF 


CELLS BETWEEN EVAPORATOR AND CONDENSER ENDS 
IF ( I . NE . 1 . AND . I . NE . I MAX ) THEN 


LIQUID/VAPOR INTERFACE CELLS 

IF (J.EQ.JV) THEN 
IF (I.LE.IE) THEN 

™ f D0/ ( TVESf1 **l • 5 >-l . 5D0/C TVESn*#2 . 5 )*TVES 
IXJ R !H = i • 5D0/DSQRT(TVESn>-0 . 5D0/ ( TVESP1XX1 . 5 )#TVES 
PGh=PGAE 
ELSE 

IwISH”? • 5D0/CTVCSM**! . 5 )-l . 5D0/( TVCSMXK2 . 5 >*TVCS 
^ T ^- 5D0/DS0R T ( T VC S” J -» ' 5D0/C TVCSPIXXl . 5 )*TVCS 

. ENDIF 
A( K, 1 )=0 . DO 

ACK, 2 )=A2-AlxPG4#TVTRm 

A(K,3)=A1+A3 

Dun=-A1XPGAXTVTRPI2 

_ )=B1 * T ( 1-1 * J >+B2*T ( 1 . J )+B3#T ( I + 1 , J )+DUf1 

ENDIF 


INTERIOR CELLS 


IF ( J . NE . J V . AND . J . NE . JU) THEN 
ACK, 1 )=A1 
ACK.2)=A2 


A(K, 3 )=A3 

B CK)=B1XT(I-1, J >+B2*T(I, J)+B3*T(I+1 
ENDIF 


,J) 
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PIPE AXIAL WALL CELLS 

IF (J.EQ.JW) THEN 
ACK, 1 )=A1+A3 

ACK,2)=A2-A3*4.D0XCA*SRADCI )*CTWSCI )**3) 

ACK.3)=0 .DO 

BCK)=B1*TCI-1,J )+B2*T C I , J )+B3*T C 1+1 , J ) 

+A3*CA#CQSTRCI >-SRADC I )* C3 . DO* C TWS C I )#*« >+TE4) ) 

ENDIF 
ENDIF 

CONDENSER END CELLS 
IF CI.EQ.IMAX) THEN 

CORNER BETWEEN L/V INTERFACE AND CONDENSER END 
IF (J.EO.JV) THEN 

TVTRni=2 . 5D0/ ( TVCSMXX1 . 5 )-l . 5D0/( TVCSM**2 . 5 >#TVCS 
TVTRf12=l . 5D0/DSQRT CTVCSPI)— 0 . 5D0/CTVCSf1*Xl . 5 >*TVCS 
ACK,1)=0.D0 

AC K . 2 )=A2-A1#PG‘+CXTVTRM1 
A ( K , 3 )=A1+A3 
Dun=-A1XPG<+CKTVTRP12 
BCK)=CB1+B3)*TCI. J)+B2KTCI,J>+DUM 
ENDIF 

CONDENSER END BETWEEN CORNERS 

IF CJ.NE. JV.AND. J.NE.JW) THEN 
AC K , 1 >=A1 
ACK,2)=A2 
AC K , 3 )=A3 

BCK)=CB1+B3)*T C IP1P11 > J )+B2*T C IflAX, J ) 

ENDIF 

CORNER BETWEEN CONDENSER END AND CONDENSER WALL 

IF C J.EQ.JW) THEN 
ACK.1)=A1+A3 

A C K . 2 ) = A2-A3K4 . D 0*CA*SRAD C I ># C TWS C I ) «*3 ) 

AC K» 3 ) = 0 . DO 

BCK)=CB1+B3)*TCI-1, J)+B2*T(I, J) 

+ +A3*C4*CQSTRCI>-SRAD(I>*C3.D0KCTWSCI )K*4)+TE4) ) 

ENDIF 
ENDIF 
190 CONTINUE 

199 CONTINUE 
RETURN 
END 
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C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
G || SUBROUTINE SPLIT2 - ENERGY EOUATION SOLUTION x 

G ~ I*I S routi NE USES PEACEMAN-RACHFORD TEMPORAL ADI TO SOLVE THE X 

C X SECOND TIME SPLIT OF THE ENERGY EOUATION J 

c ^ 

C X * 

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

c 

SUBROUTINE SPLIT2 
IMPLICIT REALX8 (A-H.O-Z) 

PARAMETER (IE=3. 

1 I C=6 , 

2 IMAX=21 , 

3 JV=1, 

4 JL=A, 

5 JW=12 , 

6 I JMAX=IMAXXJW) 

COMMON /FLOW/U(IMAX, JV: JW), VCIMAX, JV: JW) 

/F1 - 0WS/ US( I MAX, JV : JW) , VS( I MAX, JV : JW> 

COMMON /FLOWN/UN ( I MAX, JV : JW) , VN( I MAX . J V : J W ) 

SSBSS 

SSSU ^lSI/Tii?^2"Ai“TEl N<In “- JV:JU> - TUS<In4X> 

COMMON /TEMP3/TVES, TVESM, TVEP1, TVCS, TVCSM, TVCP1 
COMMON /LUD/A(IJMAX,3),BCI JMAX) ’ 1 

G °™°N /ADIAS1/ ADIA1(JW),ADIA2(JW),ADIA3(JW) 
rnMMnw ^ A ^®f^ A DIBl ( JW) , ADIB2C JW) , ADIB3( JW) 

22««2^ ^ ADI AS2/AD lAA(JW), ADIA5CJW) , ADI A6 ( JW ) 

22Dm2 N /a DIBS2/ADIBA ( JW) , ADIB5C JW) , ADIB6CJW) 

ssij phirox - phixrh - phirxh 

ogfSSS ^tSSSSgi: ■ SIGnA - Unu - UNU - Cp 

CoffioU /S? P f«!RTO«- RH0E - RH0C - DPDTE - DPDTC - PSIE - PSIC 

XWPROP/SRAD ( I MAX ) , ALPHAW, TKW 

COMMON /RAD/R(JW).RS(JW).RS1(JW),RS2(JW),RS3CJU) RSAf ILH p «f ILn 
COMMON /R/PI . DELTA. RV. RL. RW. JLM1 . JLP1. JVPl7jVP2 JVP3 RS5( JW) 

COMMON /XD/XC IMAX) , IEM1 , IEP1 , ICM1 , ICP1 , IMM1, IMM2 
22mm 2 N /XDIf1 /EL,AL»CL,PIL,ELPH,CLPH 

rnMMOM HX ' HT * HXHALF. HRHALF. KTM 

COMMON /ITER/EPSTV.MAXITV. ITRTV 
COMMON /WICK/CW 

rniimw J * PG2 * PG3 * PGAE * P <*C. PG5E . PG5C, PG6 . PG7 

COMMON /FGROUP/FO , FI , F2 » F3 , FA , F5 . F6 F7 Ffi Ti^ 

COMMON /CGROUP/C1 , C2 ! C3 . CA ' * 13 

COMMON /VGROUP/V1. V2, V3, VA, V5 

COMMON /ZGROUP/Z1 . Z2 . Z3 . ZA . Z5 . Z6 . Z7 . Z8. Z9 . Z1 0 . Z1 1 . Z12 . Z13 

CALCULATE AT FULL TIME 

DO 599 J=JV,JW 

DO 590 1=1. IMAX 
K=I+IMAXXC J-l ) 


CALCULATE ELEMENTS 

AA=ADI AA C J )-HXHALFXUC I , J) 
A5=ADIA5( J) 

A6=ADIA6(J)+HXHALFXUCI. J) 
BA=ADIBA C J )-VS( I . J >XPHIXRH 
B5=ADIB5 ( J ) 

B6=ADIB6CJ)-VSCI. J)XPHIXRH 
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LIQUID-VAPOR INTERFACE CELLS 
IF (J.EQ.JV) THEN 

CORNER BETWEEN L-V I/F AND EVAPORATOR END 

IF CI.EQ.l) THEN 
A(K, 1 )=0 .DO 
A ( K , 2 ) = A5 
AC K , 3 )=AA+A6 

B ( K >= C B5-BAXPG4E/ ( TVES*#1 . 5 ) )#T( I . J ) 

1 +CBA+B6)*TCI.J+1>+BA*PGAE/DSQRTCTVES) 

ENDIF 

CELLS BETWEEN ENDS 

IF ( I . NE . 1 .AND . I . NE. I MAX) THEN 
ACK,1)=A<* 

AC K , 2 )=A5 
ACK,3>=A6 
IF CI.LE.IE) THEN 
TEP1P=TVES 
PGA=PGAE 
ELSE 

TErtP=TVCS 

PGA=PGAC 

ENDIF 

B(K)=( B5-BAXPG4/ C TEMP**1 . 5 ) >«T C I . J V ) 

1 + ( B4+B6 )*T ( I . JV+1 )+BA*PG4/DSQRT ( TEflP ) 

ENDIF 

CORNER BETWEEN L-V I/F AND CONDENSER END 

IF (I.EQ.inAX) THEN 
ACK,1)=AA+A6 
ACK,2)=A5 
A(K,3)=0 .DO 

B ( K >= CB5-B4XPGAC/CTVCSHX1 . 5 > >*TC I . JV > 

1 + ( B4+B6 )*T ( I , JV+1 >+BA*PGAC/DSQRT ( TVCS ) 

ENDIF 
ENDIF 

CELLS BETWEEN L/V INTERFACE AND PIPE AXIAL WALL 

IF (J.NE. JV. AND. J.NE. JW) THEN 

EVAPORATOR END CELLS 

IF (I.EQ.l) THEN 
ACK, 1 )=0 . DO 
ACK.2)=A5 
ACK, 3 >=AA+A6 

BCK)=BAKTC1, J-1)+B5*TC1.J>+B6XT(1, J+l) 

ENDIF 

CONDENSER END CELLS 

IF CI.EQ.IflAX) THEN 
A(K.1>=A4+A6 
ACK,2)=A5 
ACK,3)=0.D0 

B ( K ) =B AXT C I MAX. J-l >+B5*T ( I MAX. J >+B6*T ( I MAX, J+l ) 
ENDIF 
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c 

c 

c 


c 

c 

c 

c 

c 

c 


+ 


c 

c 

c 


INTERIOR CELLS 


IF C I . NE . 1 . AND . I . NE . I MAX ) THEN 
ACK,1)=A4 


A C K, 2 ) = A5 
AC K , 3 )=A6 

BCK) = BAXT(I, j-i )+B5*TCI, J)+B6*TC 
ENDIF 
ENDIF 


I 


. J + l ) 


PIPE AXIAL WALL BOUNDARY CELL 

IF (J.EQ.JW) THEN 

EVAPORATOR CORNER CONDITION 

IF (I.EO.l) THEN 
ACK, 1 )=0 . DO 
AC K, 2 )=A5 
ACK, 3 )=AA+A6 

BCK)=CBA+B6)*TCI, J— 1 )+B5*T C I » J ) 

CKinrc _ B6XC<+X ( OSTR ( I )+SRAD( I )X ( T( I . J )XX4-TE4 ) ) 


AXIAL WALL CONDITION 


C 

c 

c 


IF ( I . NE . 1 . AND . I . NE . I MAX ) THEN 
ACK,1)=AA 
ACK, 2 )=A5 
AC K, 3 )=A6 

BCK)=CB4+B6)*TCI, J-1)+B5*T(I, J) 

1 ENDIF -W^CCWTRCn^SRADCI^CTCI.J^-TEO) 

CORNER BETWEEN AXIAL WALL AND CONDENSER END 

IF CI.EQ.IMAX) THEN 
ACK, 1 )=AA+A6 
ACK, 2 )=A5 
ACK, 3)=0 . DO 

BCK)=CBA+B6)XTCI,J-1)+B5KTCI,J) 

1 ENDIF ~ B 6XCAX COSTRC 1 )+SRAD(I )*(T( I , J )X*4-TE4 ) ) 

ENDIF 

590 CONTINUE 

599 CONTINUE 
RETURN 
END 
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c XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C X * 


SUBROUTINE LORI 

THIS SUBROUTINE USES LOWER-UPPER DECOriPOS I T I ON TO SOLVE 
TRIDIAGONAL LINEAR SYSTEMS OF EQUATIONS OF THE FORM AX=B. 
DEFINITIONS: 

A = TRIDIAGONAL MATRIX (SUPPLIED BY CALLING PROGRAM) 

B = RESULTANT VECTOR (SUPPLIED BY CALLING PROGRAM) 

X = SOLUTION VECTOR (REPLACES B ELEMENTS IN SOLUTION) 
UMAX = NUMBER OF ELEMENTS IN VECTOR B 

THE ARRAY NOTATION A(ROW 8 , BAND 8) IS USED 

BAND 8 IS ASSIGNED FROM LOWER LEFT CORNER TO UPPER RIGHT 


c x * 

C xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

c 

SUBROUTINE LORI 
IMPLICIT REALX8 (A-H.O-Z) 

PARAMETER (IE=3, 

1 IC=6, 

2 IMAX=21 , 

3 JV=1, 

A J L=A , 

5 JW=12, 

6 I JMAX=IMAXXJW) 

COMMON /LUD/A( UMAX, 3 ) , B ( UMAX) 

IJM1=I JMAX-1 


C 

C CONSTRUCT LOWER AND UPPER MATRICES FROM ELEMENTS OF A 

C 

DO 17 1 = 1. UMAX 

IF (I.GT.l) THEN 

A(I»2)=A(I,2 )— A( 1 , 1 )XA( I— 1 » 3 ) 

ENDIF 

IF (I . LT. UMAX) THEN 
A( I , 3 )=A( I . 3 )/A(I , 2 ) 

ENDIF 
17 CONTINUE 
C 

C FORWARD SUBSTITUTION IS USED TO SOLVE THE INTERMEDIATE EQUATION LY=B. 
C 

B(1)=B(1) /A (1,2) 

DO 20 1=2, UMAX 

B(I )=(B(I>-A(I,1)XB(I-1))/A(I,2) • 


20 CONTINUE 


C 

C BACK SUBSTITUTION IS USED TO SOLVE FOR X FROM UX=Y. X IS STORED IN B 
C 

DO 25 1=1. Util 

B ( UMAX— I )=B ( I JMAX-I )-A(IJMAX— 1 , 3 )*B( I JMAX-I+1 ) 

25 CONTINUE 
RETURN 
END 
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XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX* 

* SUBROUTINE PROPS - LIQUID THERMOPHYSICAL PROPERTIES 

* 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 


XXXXXXXXXXXXXXXXXXXXXXXXX 


c 

c 

c 


SUBROUTINE PROPS 
IMPLICIT REALMS (A-H.O-Z) 

PARAMETER (IE=3, 

1 IC= 6 , 

2 IMAX=21 , 

f JV=1 , 

* JL=A , 

5 JW=12, 

* I JMAX=IMAXXJW) 

COMMON /FLOWS/US ( IMAX, JV : JW) , VS( IMAX. JV : JW) 

CQMMnw ✓ci L Si U !^ /UNCInAX ’ JV: JW >»VN(IMAX.JV:JW) 

COMMON /FLOWT/UTC IMAX, JV: JW), VTCIMAX JV'JLI) 

COMMON /UBAR/UB(IMAX);SAMMACIMAX) ,JV ' JW> 

CoSK! /liriP 2 /TVl?TV 2 ^TAV W TEl NtInAX ' JV:JW, * TWSaP,A X ) 

S ; L T S/^L V nAiI 3 V !f?a"nKr VCS ' TVCSn - TVCPI 

/ADI AS 1 /ADI A1 C JW) , ADI A2 ( JW) , ADI A3 ( JW) 
common /An?!fo /ADIB1(JU) ’ AD IB2(JW).ADIB3( JW) 
rnMMHM ✓ * £ ™ i?£ /AD 1 AA C Ji W 5 * ADI A5 ( J W ) . ADI A6 ( JW) 
rnwMHM ^SJ1 S « /ADIBA C JW) * ADIB5 ( JW) , ADIB6C JW) 

SmSS ^EJKqST^^ 

/VPROPl/vnu! VNLJ.RHOvfpVAp 1 " S ^ GP,A ' UnU * UNU ‘ CP 
cofrioN RH0E> RH0C - DPDTE> DPDTC - PSIE - PS,C 

/WPROP/SRADC IMAX) , ALPHAW, TKW 
common /?t« A f XHR * HX * HT * hxhalf, hrhalf. ictm 

“SSJ ^JIK!S STV - nAXITV - ITI!TV 

C§5SoK ^g^?J;^Jf^^- P f 5 | 4 P «C. p G6. p S 7 

COMMON /CGR0UP/C1 . C2 . C3 ! CA ' ' 13 

COMMON /VGROUP/V 1 , V2 , V 3 , VA V 5 

COMMON /2GR0UP/21 . 22 , 23 . Z4 , 25 . 26 . Z7 , Z8 . 29 . Z 1 0 , ZU , 212 , 213 
UPDATE THERMOPHYSICAL PROPERTIES USING AVERAGE VAPOR TEMP 
TAV= 0 . 5D0* ( TVEP1+TVCP1 ) 

PVAP=101325 DOXDEXPd8.832DO-13113.DO/TAV 

TR=TAi>?CR DL0G<TAV)+1 * 9777 D * A *TAV> 

TAV2=TAVXTAV 

TAV3=TAVXTAVXTAV 

SIGMA= . 23402D0-1 . D-4XTAV 

■ 0578D2-5 . 1767D-2XTAV+A . 8696D-6XTAV2 
nDn^m^ 8706 " 1 ! .D0-TR)+3. A860D6X( (l . D0-TR)XX0 2) 
DPDT=PVAPX( 13113. D0/TAV2-1 . 09A8D0/TAV+1 . 9777 D— 4 ) 
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IF (TAV LE 1644. DO) THEN 

RHOL= 1 Oil 8D0-. 22054 DOKTAV-l . 9226D-5XTAV2+5 . 6371D-9XT AV3 
DRODT=-. 22054D0-3.8452D-5XTAV+1 . 691 1D-8#T AV2 

DHDT=1 . 53139D3- . 613452D0#TAV+3 . 35513D-4XTAV2+5 . 40593D6/TAV2 
ELSE 

RHOL=RHOCRX ( 1 . DO+2 .3709D0XCC1. DO-TR )xx . 31645 ) 

1 +2.8467D-7X((TCR-TAV)*X2)) 

DRODT=-RHOCRX( . 75027D0X ( C 1 . DO-TR )*x . 68355 )/TCR 
1 +5.6934D-7* ( TCR-TAV ) ) 

DHDT=805 .8D0+304 .21D0XC ( 1 . D0/( 1 . DO-TR ) )**0 • 67773 ) 

ENDIF 

RHOV=l ,D0/CHFG/CTAV*DPDT)+1 . DO/RHOL) 

Uf1U=l . 1259D-5XCRHOLXX ( 1 . DO/3 . DO ) )*DEXP( . 74908XRHOL/TAV) 
UNU=UMU/RHOL 

DRODP=1^77587D-10/a.2773035D0-1.8267670XTR+.54946350D0XTRXTR) 

1 - ( 2 15463D-21XTAV+1 .2204 08D-24XTAV2 )XPVAP 

2 +1 . 671245D-29XTAVXPVAPXPVAP 

CP=DHDT+ ( ( DPDTXDRODP-DRODT >*TAV/RHOl-l - DO )XDPDT/RHOL 
ALPHAL=TKL/ (RHOLXCP ) 

UPDATE EVAPORATOR REGION VAPOR PROPERTIES 

PVE=1 01325 . D0XDEXPC18.832D0-13113.D0/TVEP1 
1 -1 . 0948D0XDLOGCTVEP1 ) + l . 9777D— 4XTVEP1 ) . 

DPDTE=PVEX (13113. D0/CTVEP1XTVEP1 )-l - 0948D0/TVEP1+1 . 9777D-4 ) 
RHOE= 1 . DO/ ( HFG/ C TVEP1XDPDTE ) + l . DO/RHOL ) 

UPDATE CONDENSER REGION VAPOR PROPERTIES 

PVC=1 01325 . DOXDEXPC 18 . 832D0-131 13 . D0/TVCP1 
1 -1 0948D0XDLOGCTVCP1 ) + l . 9777D— 4XTVCP1 > , . 

DPDTC=PVCX C 13113 . D0/CTVCP1XTVCP1 )-l . 0948D0/TVCP1+1 . 9777D-4 ) 
RHOC= 1 . DO/ C HFG/ C TVCP1XDPDTC ) + 1 . DO/RHOL ) 

UPDATE PROPERTY DEPENDENT PARAMETER GROUPS 


PG1=V3XHFG 

PG2=V4XHFG 

PG3=RHOLXHFGXCUI 

PG4E=C3XHFGXHFGXRH0E/TKL 

PG4C=C3XHFGXHFGXRH0C/TKL 

PG5E=C2XHFGXRH0E/RH0L 

PG5C=C2XHFGXRH0C/RH0L 

PG6=V5XHFG 

PG7=RHOLXRHOLX C 1 . DO/RHOV-1 . DO/RHOL ) 

RETURN 

END 
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XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
* SUBROUTINE REGIM1 - UNIFORM VAPOR CONDITION X 
X * 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


c 

c 

c 


200 

c 

c 

c 


SUBROUTINE REGIM1 (ISPLIT) 

IMPLICIT REALX8 (A-H.O-Z) 

PARAMETER CIE=3, 

1 IC=6, 

2 IMAX=21 , 

3 JV=1, 

A JL=4, 

5 JW=12, 

6 I JMAX= I MAXX JUI ) 

COMMON /FLOW/UCIMAX, JV: JW), VCIMAX.JV: JW) 

COMMON /FLOWS/US ( I MAX, JV : JW) , VS( IMAX, JV: JW) 

COMMON /FLOWN/UN ( IMAX, JV: JW) , VN(IMAX, JV: JW) 

COMMON /FLOWT/UT (IMAX, J V : JW ) , VT ( IMAX, JV : JW) 

COMMON /UB AR/UB ( I MAX ) , GAMMA ( I MAX ) 

COMMON /TEMP1/TCIMAX, JV: JW) , TN( IMAX. JV: JW) , TWSC IMAX) 

COMMON /TEMP2/TV1, TV2, TAV. TEA iwiUHAXJ 

COMMON /TEMP3/TVES, TVESM, TVEP1 , TVCS, TVCSM, TVCP1 
COMMON /LUD/ AC UMAX, 3 ) , B ( UMAX) 

COMMON /AD IASI /ADI A1 ( JW) , ADI A2( JW) , ADI A3( JW) 

COMMON /ADIBS1/ADIB1 (JW) , ADIB2( JW). ADIB3( JW) 

rnwwnM n^ /ADI A4 C JW 5 ' ADI A5 ( JW) * ADI A6 ( JW > 

COMMON /ADIBS2/ ADIB4CJW), ADIB5C JW) , ADIB6C JW) 

c§ss^f;^^;jnK;r=ii3 IX0R - PHiR0x ' PHixRH - pHiRxH 

SoSS ./tpS p ^!Si: SSi: SIGm - unu - UNU - cp 

comm /CRI t^tcrUhocr ■ RH ° E ' ■ RHDe - : DPDTE - DPDTC - PSIE • PSI c 

COMMON /WPROP/SRAD (IMAX ) , ALPHAW, TKW 

COMMON /RAD/RC JW) , RS( JW) , RSI ( JW) , RS2( JW) , RS3 ( JW) RSAfJUl p« f 
COMMON /R/PI , DELTA, RV, RL , RW. JLM1 JLP1 , JVP1 .j5p2. JVp5 

SSS / ^il? AX J‘l em ‘ 1EP1 - lcni ‘ icpi,inrii,iMM2 JVP3 
COMMON /XDIM/EL . AL , CL , PI L , ELPH.CLPH 

rnMMnlll HX ' HT * HXHALF, HRHALF. KTP1 

COMMON /I TER/EPSTV , MAXI TV , I TRT V 
COMMON /WICK/CW 

rnMMnM * t G2 ’ PG3 ' PG4E * PG4C * PG5E. PG5C, PG6 , PG7 

COMMON /FGROUP/F 0 , FI , F2, F3, FA , F5,F6,F7,F8, T13 
COMMON /CGROUP/C1 , C2 , C3 , CA 
COMMON /VGROUP/V1, V2.V3, VA , V5 

COMMON /2GR0UP/Z1 , Z2 » Z3 , 2A , Z5 , 26 , 27 , Z8 . Z9 , 21 0 , Z1 1 , 212 , Z13 
SUM SURFACE TEMPERATURE 


SUM=0 . 5D0X (T( 1 , JV)+T(IMAX, JV) ) 
DO 200 1=2, IMM1 
SUM=SUM+T (I , JV) 

CONTINUE 


SET PARAMETERS FOR NEWTONS METHOD 

IF (ISPLIT. EQ. 1) TV1=TVES 
IF (ISPLIT.EQ.2) TV1=TVEP1 
BETA=HXXSUM 
DRATI 0=-RHOE/RHOL 
TERMN=-RHOE*HFG/DPDTE 
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USE NEWTONS METHOD TO FIND VAPOR TEMPERATURE 

DO 600 M=l, MAXITV 

TVM1H=1 .D0/DSQRTCTV1) 

TVM3H=1.D0/(TV1*«1.5) 

AETA=1.D0/(1. D0+PG6#TVM3H* ( TV1*PI L-BETA ) ) 

FTV=TV1# ( AETA+DRATI 0)+TERMN 

FPRIME= AETA# ( 1 . D0-AETA*PG6*0 . 5D0# (3 . D0*BETA*TVM3H- 
1 PIL#TVM1H) )+DRATIO 

DELTV=-FTV/FPRIME 
TVl=TVl+0 . 8D0#DELTV 

itrtv=pi 

IF CDABSCDELTV) .LT.EPSTV) GOTO 700 
600 CONTINUE 

700 IF (M.EQ. MAXITV) WRITE (6,9050) M,KTM 
TV2=TV1 

9050 FORMAT ( /////l OX, * ***** CAUTION *#**#’ 

1/10X, ’VAPOR TEMPERATURE DID NOT CONVERGE 1 
2/15X, ’ ITERATIONS = ’,13 
3/16X, ’TIME STEP = ’,13) 

RETURN 

END 
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c 
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SUBROUTINE REGIP12( ISPLIT) 

IMPLICIT REAl#8 (A-H.O-Z) 

PARAMETER (IE=3, 

1 IC=6, 

2 IMAX=21 , 

3 JV=1. 

* Jl=4, 

f JW=12, 

6 I JMAX=IMAX#JW) 

COMMON /FLOW/UC IMAX, JV : JW) V C IMAX JV* Itn 
COMMON /F L OWS/US ( I MAX . J V : J W ) . VS (I MAX ’ jS? m > 

COMMON /FL OWN/UN ( IMAX, JV : JW) VN( IMAX* JV • JLM 

ass 

COMMON /TB?2^V1 ?TV2^TAV W TeI N C JV = JW) * TWS ( InAX J 

SK^^^fSiiSjTV'S.TVCSn.TVCP! 

SBSB 

SBS 

p ag@aaEc==r™ - 

S Sis 

Hiassapcr-" 

COMMON /FGROUP/FO ’ F1 2 F2^F3^FA^F5**FA ’p? 5 ^!*^ 50 ' PG6 * PG7 
COMMON /CGROUP/C1 fclfclfcl ’ ’ F6 ' F7 ‘™‘ 113 

COMMON /VGROUP/V1, V2.V3. V4 V5 

COnnON /2GR0UP/21 ■ 22, 23, 2^ izS, 26. 27. 28, 29, 210, Zll, 212. 213 


FIND EVAPORATOR VAPOR TEMPERATURE 


IF (ISPLIT. E0.1 
IF (ISPLIT. EQ. 2 
SUM=T(1, JV)#0.51 
DO 200 1=2, IE 


TV1=TVES 

TV1=TVEP1 


SUM=SUM+T ( I , JV) 
CONTINUE 
BETA=HX#SUM 


RATIO=-RHOE/RHOL 
TERMl=-PSIE/RHOE 
TERM2=— RHOEKHFG/DPDTE 
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DO 300 M=1,MAXITV 

TM1H= 1 . DO/DSQRT ( TV1 ) 

TP13H=1 . D0/TV1KK1 . 5 

AETA=1 . D0/( 1 . D0-PG1*ELPHXTM1H+PG1XBETA*TM3H+TERM1 ) 
FTV=AETA*TV1+RATI0*TV1+TERM2 

FPRin=AETA* C 1 . DO-O . 5D0*AETAXPGl*(ELPH*TMlH-3 . D0#BETA*TM3H ) ) 
+ +RATIO 

DELTV=-FTV/FPRIM 
TV1 = TV1 + 0 . 8D0KDELTV 
ITRTV=n 

IF (DABS(DELTV) . LT . EPSTV) GOTO A00 
300 CONTINUE 

IF Cn.EQ.nAXITV) -WRITE (6,9050) M, KTM 
FIND CONDENSER VAPOR TEMPERATURE 


A00 IF (ISPLIT.EQ. 1) TV2=TVCS 
IF (ISPLIT.EQ. 2) TV2=TVCP1 
SUM=T(IMAX, JV)X0.5D0 
DO 500 I=IE+1 , IMAX-1 
SUM=SUM+T(I, JV) 

500 CONTINUE 

BETA=HX#SUM 
RATIO=-RHOC/RHOL 
TERMl=PSIC/RHOC 
TERM2=-RH0C*HFG/DPDTC 
DO 600 M= 1 , MAXI TV 

TM1H= 1 . D0/DSQRT ( TV2 ) 

TM3H=1 . D0/TV2XX1 .5 

AETA*1 D0/( 1 . D0-PG2#CLPH*TM1H+PG2XBETA*TM3H+TERM1 ) 
FTV=AETAXTV2+RATIO*TV2+TERM2 

FPRIM=AETA#( 1 . DO-O . 5D0XAETAXPG2XCCLPHXTM1H-3 . D0XBETAXTM3H) ) 
+ +RATIO 

DELTV=-FTV/FPRIM 
TV2=TV2+0 . 8D0XDELTV 


ITRTV=M 

IF (DABS(DELTV).LT. EPSTV) GOTO 700 
600 CONTINUE 

700 IF (M.EQ. MAXITV) WRITE (6,9050) M.ICTM 
RETURN 

9050 FORMAT C/////10X, • «***« CAUTION *«*««• 
1/1 OX, ’VAPOR TEMPERATURE DID NOT CONVERGE 


2/15X, 'ITERATIONS = ’.13 
3/16X, ’TIME STEP * ’.13) 
END 
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C KXKKKKXK***##****#*#*****##*****#****#****##******###*#****^ 

<r * SUBROUTINE SPEED - LIQUID VELOCITIES * 

L * 

c x * 

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXHXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXX 


SUBROUTINE SPEED( ISPLIT) 

IMPLICIT REALX8 (A-H.O-Z) 

PARAMETER (IE=3, 

1 I C=6 , 

2 IMAX=21 , 

3 JV=1, 

A JL=A, 

5 JW=12, 

6 IJMAX=IMAXXJW) 

COMMON /F L 0 W/U C I MAX , J V : J W ) , V C I MAX » J V : J W ) 

COMMON /FLOWS/US ( I MAX, JV : JW) , VSC IMAX, JV : JU) 

COMMON /FLOWN/UN (IMAX, JV : JW ) » VN ( IMAX. JV : JW) 

rnMMHM /up ?d / u» T ( IMAX ' JV ; JW) , VTC IMAX, JV : JW) 

COMMON /UB AR/UB ( I MAX ) , GAMMA ( I MAX ) 

cSKS ^u^AajSi!??? B n ; IK; Tvcs - Tvcs "- TVCP1 

/ADIAS1/AD IA1(JW),ADIA2(JW).ADIA3(JW) 

< A DIBS1/ADIB1CJW),ADIB2CJW).ADIB3CJW) 

rnMMnw /An^li/A^ AA(JW) * ADIA5CJW),ADIA6(JW) 

COMMON /ADIBS2/ADIB4 ( JW) . ADIB5( JW) , ADIB6 ( JW) 

SSK ;^?^r!K!^IU ,xor - phirox - ,,hi><rh ’ phirxh 

SSS x t p RS p ^nt : ™ • si ««- . cr 

comm ^ru P ?5I1hocr ' RH0E - RH0C - DPI>TE - DPDTC - P£IE - PSIC 

COMMON /WPROP/SR AD C I MAX ) , ALPHAW, TKW 

COriMON /R/Pl J p, ) ’ dS 1 ttK? • Rp£ ( JW ) , RS3 C JW) . RSA ( J W ) . RS5 ( JW) 

P^TA, RV,RL, RW, JLttl * JLP1 , JVP1 , JVP 2 , JVP3 
COnMON /XD/XCIMAX), I EMI, IEP 1 , ICM 1 , ICP 1 , IMM 1 , IMM 2 
COMMON /XDIM/EL, AL, CL,PIL,ELPH,CLPH 

/?^I^I» H 5* HX * «T. HXHALF. HRHALF, KTM 
COMMON /I TER/EPSTV, MAXI TV, ITRTV 
COMMON /WICK/CW 

COMMON /PROPG/PG1 , PG2, PG3 , PGAE, PG4C, PG5E, PG5C PGA PG7 
COMMON /FGROUP/FO ,F1,F2,F3, Fa7f5 . F6 F7.F87T13 * ’ 

COMMON /CGROUP/C1 , C2 , C3, C4 
COMMON /VGROUP/V1, V2.V3, VA, V5 

COMMON /2GR0UP/21 . Z2 . Z3 . ZA , Z5 , 26 . Z7 . 28 . Z9 , Z1 0 , Z1 1 . Zl 2 . Z1 3 
FIND FIRST TIME SPLIT L-V I/F RADIAL VELOCITIES 


100 

150 


IF 


( ISPLIT . EQ . 1 ) THEN 
DO 100 1=1. IE 

VTC I , JV)=(TVES— T(I 
CONTINUE 


. JV) )XPG5E/TVES*X1 


DO 150 I=IEP1 , IMAX 

C0N^aE JV)=CTVCS ' Ta ' JV)>XPG5C/TVCS,<K1 

END IF 


.5 

.5 
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c 

C FIND SECOND TIME SPLIT L-V I/F RADIAL VELOCITIES 

C 

IF (ISPLIT.EQ.2) THEN 
DO 200 1=1, IE 

VTCI, JV)=CTVEP1-TCI , JV) )*PG5E/TVEP1 *#1 .5 
200 CONTINUE 

DO 250 I = I E P 1 , I MAX 

VT(I , JV)=CTVCP1-T(I , JV) )*PG5C/TVCP1**1 .5 
250 CONTINUE 

ENDIF 


C 

C FIND MEAN AXIAL VELOCITIES 

C 

SUM=VT ( 1 , JV ) 

DO AO 0 1=2, IMM1 

UBCI )=C1*(SUM+VT(I, JV) ) 

SUM=SUM+2 . DQ#VT C I , J V ) 

AO 0 CONTINUE 
C 

C FIND GRID VELOCITIES 

C 

DO 700 1=2, IMM1 

IF (UBCI).NE.O.DO) THEN 

IF (VTCI, JVKNE.O.DO) THEN 

DU2DX= ( UB ( I -1 )-2 . D0#UB ( I )+UB ( 1+1 ) )/CHX*HX ) 

S1=F7*VT(I, JV5/UNU 

S2=FA+F5#DU2DX/UB ( I )+F6#VT ( I , J V )/UNU 
S3=F1+F2*DU2DX/UB ( I )+F3#VT ( I , J V )/UNU 
DD=(-S2-DSQRT (S2#S2-A . D0*S1*S3) )#0.5D0/S1 
ELSE 

DD=0 . DO 
ENDIF 

BB=6.D0+T13*DD 

CC=-BB-DD 

ENDIF 

DO 600 J=JVPl,JLf11 

UT ( I , J )=UB C I )#(BB*RS ( J >+2 . D0*CC*RS2 ( J )+3 . D0*DD#RS3 C J ) ) 

VTCI , J) =CU*VT (I , JV)K( 1 . D0-F8* ( BBXF0KRS2 ( J )+ ( BB+CC*F0 ) 

1 *RS3 ( J )+ ( CC+DDKFO >*RSA C J >+DD*RS5 ( J ) ) >#RV/R ( J ) 

IF (ISPLIT.EQ.2) THEN 

GAMMA ( I )=UMU* ( Z1 *BB+Z2*CC+Z3*DD )+RHOL#VT ( I , J V ) 

1 X ( ZA*BBHBB+Z5#BB*CC+Z6#BB*DD+Z7*CC#CC 

2 +Z8XCCXDD+Z9XDDXDD ) 

ENDIF 

600 CONTINUE 

700 CONTINUE 
RETURN 
END 
C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C INPUT DATA 

c Nun 

C NCASE HX HR RV 

C PR UR HT TE 

C TV 

C 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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FIND SECOND TIME SPLIT L-V I/F RADIAL VELOCITIES 
IF CISPLIT.EQ.2) THEN 

D ° VT(l! JV j=(TVEPl-T(I, JV) )*PG5E/TVEP1**1.5 
200 CONTINUE 

DO 250 I=IEP1 > I MAX . c 

VTCI.JV)=<TVCP1-TCI.JV)>*PG5C/TVCP1**1.5 
250 CONTINUE 

ENDIF 


C 

c 

c 


400 

C 

c 

c 


600 

700 


FIND MEAN AXIAL VELOCITIES 

SUM=VT ( 1 . JV ) 

DO 400 1 = 2 > IMM1 

UB(I )=C1* ( SUM+VT ( I , JV ) ) 
SUM=SUM+2 . D0XVT ( I , JV ) 
CONTINUE 

FIND GRID VELOCITIES 


DO 700 1=2 , IMM1 

IF (UB( I ) . NE. 0 .DO ) THEN 

IF ( VT( I . JV) .NE. 0 .DO ) THEN 

DU2DX=(UB C I — 1 )-2 . D0*UB(I )+UB(I + l ) )/(HX*HX) 


S1=F7*VT(I, JV)/UNU 
S2=F4+F5XDU2DX/UB ( I 
S3=F1+F2*DU2DX/UB(I 
DD=C-S2-DSQRT(S2*S2 


)+F6*VT(I.JV)/UNU 
)+F3«VT<I. JV)/UNU 
-4 .D0*S1*S3> )*0 .5D0/S1 


ELSE 

DD=0 . DO 
ENDIF 

BB=6 . D0+T13XDD 
CC=-BB-DD 


1 

2 


ENDIF 

D ° UT(I^J)=UBCI >*CBBXRS(J>+2.D0XCCXRS2(J>+3.D0*DDKRS3( J 

VT ( I . J ) =CW*VT C I . J V )X ( 1 . DO-F8* ( BB*F0*RS2 C J )+ ( BB+CC*F0 

XRS3C J )+(CC+DD#F0)XRS4(J>+DD*RS5( J > 5 )*RV/R( J 

1 ^ GAMMA ( I j =UMU* ( Z1*BB+Z2*CC+Z3XDD )+RHOL#VT C I , JV ) 

X(Z4XBB*BB+Z5XBB#CC+Z6XBBXDD+Z7XCC*CC 

+Z8*CC*DD+Z9XDDXDD) 


ENDIF 

CONTINUE 

CONTINUE 

RETURN 

END 


)) 

) 

) 




c 

c 

c 

c 

c 


INPUT DATA 
NUfl 


NCASE 

PR 

TV 


HX 

UR 


HR 

HT 


RV 

TE 
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